home *** CD-ROM | disk | FTP | other *** search
/ Complete Linux / Complete Linux.iso / xwindows / demos / xfract_1.z / xfract_1 / xfractint-1.06 / calcfrac.c < prev    next >
C/C++ Source or Header  |  1992-09-28  |  117KB  |  4,191 lines
  1. /*
  2. CALCFRAC.C contains the high level ("engine") code for calculating the
  3. fractal images (well, SOMEBODY had to do it!).
  4. Original author Tim Wegner, but just about ALL the authors have contributed
  5. SOME code to this routine at one time or another, or contributed to one of
  6. the many massive restructurings.
  7. This module is linked as an overlay, use ENTER_OVLY and EXIT_OVLY.
  8. The following modules work very closely with CALCFRAC.C:
  9.   FRACTALS.C    the fractal-specific code for escape-time fractals.
  10.   FRACSUBR.C    assorted subroutines belonging mainly to calcfrac.
  11.   CALCMAND.ASM    fast Mandelbrot/Julia integer implementation
  12. Additional fractal-specific modules are also invoked from CALCFRAC:
  13.   LORENZ.C    engine level and fractal specific code for attractors.
  14.   JB.C        julibrot logic
  15.   PARSER.C    formula fractals
  16.   and more
  17.  -------------------------------------------------------------------- */
  18.  
  19. #include <stdio.h>
  20. #include <string.h>
  21. #include <stdlib.h>
  22. #include <float.h>
  23. #ifndef XFRACT
  24. #include <dos.h>
  25. #endif
  26. #include <limits.h>
  27. #include "fractint.h"
  28. #include "fractype.h"
  29. #include "mpmath.h"
  30. #include "targa_lc.h"
  31. #include "prototyp.h"
  32.  
  33. /* routines in this module    */
  34.  
  35. static void perform_worklist(void);
  36. static int  OneOrTwoPass(void);
  37. static int  _fastcall StandardCalc(int);
  38. static int  _fastcall potential(double,int);
  39. static void decomposition(void);
  40. static int  bound_trace_main(void);
  41. static void step_col_row(void);
  42. static int  _fastcall boundary_trace(int,int);
  43. static int  _fastcall calc_xy(int,int);
  44. static int  _fastcall fillseg1(int,int,int,int);
  45. static int  _fastcall fillseg(int,int,int,int);
  46. static void _fastcall reverse_string(BYTE *,BYTE *,int);
  47. static int  solidguess(void);
  48. static int  _fastcall guessrow(int,int,int);
  49. static void _fastcall plotblock(int,int,int,int);
  50. static void _fastcall setsymmetry(int,int);
  51. static int  _fastcall xsym_split(int,int);
  52. static int  _fastcall ysym_split(int,int);
  53. static void set_Plasma_palette(void);
  54. static int _fastcall adjust(int xa,int ya,int x,int y,int xb,int yb);
  55. static void _fastcall subDivide(int x1,int y1,int x2,int y2);
  56. static int _fastcall new_adjust(int xa,int ya,int x,int y,int xb,int yb);
  57. static int _fastcall new_subD (int x1,int y1,int x2,int y2, int recur);
  58. static void verhulst(void);
  59. static void Bif_Period_Init(void);
  60. static int  _fastcall Bif_Periodic(int);
  61. static void set_Cellular_palette(void);
  62.  
  63. /**CJLT new function prototypes: */
  64. static int tesseral(void);
  65. static int _fastcall tesschkcol(int,int,int);
  66. static int _fastcall tesschkrow(int,int,int);
  67. static int _fastcall tesscol(int,int,int);
  68. static int _fastcall tessrow(int,int,int);
  69.  
  70. extern struct complex initorbit;
  71. extern char useinitorbit;
  72. struct lcomplex linitorbit;
  73. extern int showdot;
  74. extern unsigned int decoderline[];
  75. extern int overflow;
  76. long lmagnitud, llimit, llimit2, lclosenuff, l16triglim;
  77. struct complex init,tmp,old,new,saved;
  78. extern int biomorph,usr_biomorph;
  79. extern struct lcomplex linit;
  80. extern int basin;
  81. extern int cpu;
  82. extern char savename[80];   /* save files using this name */
  83. extern int resave_flag;
  84. extern int started_resaves;
  85. extern int dotmode;
  86.  
  87. int color, oldcolor, row, col, passes;
  88. int realcolor;
  89. int iterations, invert;
  90. double f_radius,f_xcenter, f_ycenter; /* for inversion */
  91. extern double far *dx0, far *dy0;
  92. extern double far *dx1, far *dy1;
  93.  
  94. extern int LogFlag;
  95. extern BYTE far *LogTable;
  96. extern int rangeslen;
  97. extern int far *ranges;
  98.  
  99. void (_fastcall *plot)(int,int,int) = putcolor;
  100.  
  101. extern double inversion[];        /* inversion radius, f_xcenter, f_ycenter */
  102. extern int    xdots, ydots;        /* coordinates of dots on the screen  */
  103. extern int    sxdots,sydots;
  104. extern int    sxoffs,syoffs;
  105. extern int    colors;         /* maximum colors available */
  106. extern int    andcolor;        /* colors-1         */
  107. extern int    inside;         /* "inside" color to use    */
  108. extern int    fillcolor;         /* "fillcolor" color to use    */
  109. extern int    outside;        /* "outside" color to use   */
  110. extern int    finattract;
  111. double        min_orbit;        /* orbit value closest to origin */
  112. int        min_index;        /* iteration of min_orbit */
  113. extern int    maxit;            /* try this many iterations */
  114. extern int    fractype;        /* fractal type */
  115. extern char    stdcalcmode;        /* '1', '2', 'g', 'b' */
  116. extern int    debugflag;        /* for debugging purposes */
  117. extern    int    diskvideo;        /* for disk-video klooges */
  118. extern int    calc_status;        /* status of calculations */
  119. extern long    calctime;        /* total calc time for image */
  120.  
  121. extern int rflag, rseed;
  122. extern int decomp[];
  123. extern int distest,distestwidth;
  124.  
  125. extern double    param[];        /* parameters */
  126. extern int    potflag;        /* potential enabled? */
  127. extern double    potparam[];        /* potential parameters */
  128. extern int    pot16bit;        /* store 16bit continuous potential */
  129. extern long    far *lx0, far *ly0; /* X, Y points */
  130. extern long    far *lx1, far *ly1; /* X, Y points */
  131. extern long    fudge;            /* fudge factor (2**n) */
  132. extern int    bitshift;        /* bit shift for fudge */
  133. extern long    delmin;         /* for periodicity checking */
  134.  
  135. extern double xxmin,xxmax,yymin,yymax,xx3rd,yy3rd; /* corners */
  136. extern long   xmin, xmax, ymin, ymax;           /* integer equivs */
  137. extern long   delx,dely;               /* X, Y increments */
  138. extern double delxx,delxx2,delyy,delyy2;
  139. double deltaX, deltaY;
  140. double magnitude, rqlim, rqlim2, rqlim_save;
  141. extern struct complex parm,parm2;
  142. int (*calctype)();
  143. double closenuff;
  144. int pixelpi; /* value of pi in pixels */
  145. unsigned long lm;        /* magnitude limit (CALCMAND) */
  146. extern long linitx,linity;    /* in calcmand */
  147. extern unsigned long savedmask; /* in calcmand */
  148. extern int TranspSymmetry;
  149.  
  150. /* ORBIT variables */
  151. int    show_orbit;            /* flag to turn on and off */
  152. int    orbit_ptr;            /* pointer into save_orbit array */
  153. int far *save_orbit;            /* array to save orbit values */
  154. int    orbit_color=15;         /* XOR color */
  155.  
  156. int    ixstart, ixstop, iystart, iystop;    /* start, stop here */
  157. int    symmetry;       /* symmetry flag */
  158. int    reset_periodicity; /* nonzero if escape time pixel rtn to reset */
  159. int    kbdcount, max_kbdcount;    /* avoids checking keyboard too often */
  160.  
  161. extern    int    integerfractal;     /* TRUE if fractal uses integer math */
  162.  
  163. char far *resume_info = NULL;        /* pointer to resume info if allocated */
  164. int resuming;                /* nonzero if resuming after interrupt */
  165. int num_worklist;            /* resume worklist for standard engine */
  166. struct workliststuff worklist[MAXCALCWORK];
  167. int xxstart,xxstop;            /* these are same as worklist, */
  168. int yystart,yystop,yybegin;        /* declared as separate items  */
  169. int workpass,worksym;            /* for the sake of calcmand    */
  170.  
  171. extern long timer_interval;        /* timer(...) total */
  172.  
  173. VOIDFARPTR typespecific_workarea = NULL;
  174.  
  175. static double dem_delta, dem_width;    /* distance estimator variables */
  176. static double dem_toobig;
  177. #define DEM_BAILOUT 535.5  /* (pb: not sure if this is special or arbitrary) */
  178.  
  179. /* variables which must be visible for tab_display */
  180. int got_status; /* -1 if not, 0 for 1or2pass, 1 for ssg, 2 for btm, 3 for 3d */
  181. int curpass,totpasses;
  182. int currow,curcol;
  183.  
  184. /* static vars for solidguess & its subroutines */
  185. static int maxblock,halfblock;
  186. static int guessplot;            /* paint 1st pass row at a time?   */
  187. static int right_guess,bottom_guess;
  188. #define maxyblk 7    /* maxxblk*maxyblk*2 <= 4096, the size of "prefix" */
  189. #define maxxblk 202  /* each maxnblk is oversize by 2 for a "border" */
  190.              /* maxxblk defn must match fracsubr.c */
  191. /* next has a skip bit for each maxblock unit;
  192.    1st pass sets bit  [1]... off only if block's contents guessed;
  193.    at end of 1st pass [0]... bits are set if any surrounding block not guessed;
  194.    bits are numbered [..][y/16+1][x+1]&(1<<(y&15)) */
  195. extern unsigned int prefix[2][maxyblk][maxxblk]; /* common temp */
  196. /* size of next puts a limit of MAXPIXELS pixels across on solid guessing logic */
  197. #ifndef XFRACT
  198. extern char dstack[4096];               /* common temp, two put_line calls */
  199. extern char suffix[4096];             /* cellular, put_line/get_line calls */
  200. extern unsigned int prefix[2][maxyblk][maxxblk]; /* common temp */
  201. #else
  202. BYTE dstack[4096];              /* common temp, two put_line calls */
  203. unsigned int prefix[2][maxyblk][maxxblk]; /* common temp */
  204. #endif
  205.  
  206. int nxtscreenflag; /* for cellular next screen generation */
  207. int    attractors;            /* number of finite attractors  */
  208. struct complex    attr[N_ATTR];        /* finite attractor vals (f.p)  */
  209. struct lcomplex lattr[N_ATTR];        /* finite attractor vals (int)  */
  210. int    attrperiod[N_ATTR];        /* period of the finite attractor */
  211.  
  212. /***** vars for new btm *****/
  213. enum direction {North,East,South,West};
  214. enum direction going_to;
  215. int trail_row, trail_col;
  216.  
  217. #ifndef sqr
  218. #define sqr(x) ((x)*(x))
  219. #endif
  220.  
  221. #ifndef lsqr
  222. #define lsqr(x) (multiply((x),(x),bitshift))
  223. #endif
  224.  
  225. /* -------------------------------------------------------------------- */
  226. /*        These variables are external for speed's sake only      */
  227. /* -------------------------------------------------------------------- */
  228. extern struct lcomplex lold,lnew,lparm,lparm2;     /* added "lold" */
  229.  
  230. int periodicitycheck;
  231. extern long ltempsqrx,ltempsqry;
  232. extern double tempsqrx,tempsqry;
  233. extern LCMPLX ltmp;
  234. extern int display3d;
  235.  
  236.  
  237. void calcfrac_overlay() { }    /* for restore_active_ovly */
  238. extern int first_err;   /* flag for math errors */
  239.  
  240. /******* calcfract - the top level routine for generating an image *******/
  241.  
  242. int calcfract()
  243. {
  244.    ENTER_OVLY(OVLY_CALCFRAC);
  245.  
  246.    first_err = 1;
  247.    attractors = 0;        /* default to no known finite attractors  */
  248.    display3d = 0;
  249.    basin = 0;
  250.  
  251.    init_misc();  /* set up some variables in parser.c */
  252.  
  253.    /* following delta values useful only for types with rotation disabled */
  254.    /* currently used only by bifurcation */
  255.    if (integerfractal)
  256.    {
  257.       distest = 0;
  258.       deltaX = (double)lx0[     1] / fudge - xxmin;
  259.       deltaY = yymax - (double)ly0[     1] / fudge;
  260.    }
  261.    else
  262.    {
  263.       deltaX = dx0[     1] - xxmin;
  264.       deltaY = yymax - dy0[     1];
  265.    }
  266.  
  267.    parm.x   = param[0];
  268.    parm.y   = param[1];
  269.    parm2.x  = param[2];
  270.    parm2.y  = param[3];
  271.  
  272.    if (LogFlag && colors < 16) {
  273.       static char far msg[]={"Need at least 16 colors to use logmap"};
  274.       stopmsg(0,msg);
  275.       LogFlag = 0;
  276.       }
  277.    if (LogFlag || rangeslen)
  278.       if (!(LogTable = farmemalloc((long)maxit + 1))) {
  279.       static char far msg[]={"Insufficient memory for logmap/ranges with this maxiter"};
  280.       stopmsg(0,msg);
  281.       }
  282.       else if (rangeslen) {
  283.      int i,k,l,m,numval,flip,altern;
  284.      i = k = l = 0;
  285.      while (i < rangeslen) {
  286.         m = flip = 0;
  287.         altern = 32767;
  288.         if ((numval = ranges[i++]) < 0) {
  289.            altern = ranges[i++];    /* sub-range iterations */
  290.            numval = ranges[i++];
  291.            }
  292.         if (numval > maxit || i >= rangeslen)
  293.            numval = maxit;
  294.         while (l <= numval)  {
  295.            LogTable[l++] = k + flip;
  296.            if (++m >= altern) {
  297.           flip ^= 1;        /* Alternate colors */
  298.           m = 0;
  299.           }
  300.            }
  301.         ++k;
  302.         if (altern != 32767) ++k;
  303.         }
  304.      }
  305.       else
  306.      SetupLogTable();
  307.  
  308.    lm = 4L << bitshift;         /* CALCMAND magnitude limit */
  309.  
  310.    /* ORBIT stuff */
  311.    save_orbit = (int far *)((double huge *)dx0 + 4*MAXPIXELS);
  312.    show_orbit = 0;
  313.    orbit_ptr = 0;
  314.    orbit_color = 15;
  315.    if(colors < 16)
  316.       orbit_color = 1;
  317.  
  318.    if(inversion[0] != 0.0)
  319.    {
  320.       f_radius      = inversion[0];
  321.       f_xcenter   = inversion[1];
  322.       f_ycenter   = inversion[2];
  323.  
  324.       if (inversion[0] == AUTOINVERT)  /*  auto calc radius 1/6 screen */
  325.      inversion[0] = f_radius = min(fabs(xxmax - xxmin),
  326.          fabs(yymax - yymin)) / 6.0;
  327.  
  328.       if (invert < 2 || inversion[1] == AUTOINVERT)  /* xcenter not already set */
  329.       {
  330.      inversion[1] = f_xcenter = (xxmin + xxmax) / 2.0;
  331.      if (fabs(f_xcenter) < fabs(xxmax-xxmin) / 100)
  332.         inversion[1] = f_xcenter = 0.0;
  333.       }
  334.  
  335.       if (invert < 3 || inversion[2] == AUTOINVERT)  /* ycenter not already set */
  336.       {
  337.      inversion[2] = f_ycenter = (yymin + yymax) / 2.0;
  338.      if (fabs(f_ycenter) < fabs(yymax-yymin) / 100)
  339.         inversion[2] = f_ycenter = 0.0;
  340.       }
  341.  
  342.       invert = 3; /* so values will not be changed if we come back */
  343.    }
  344.  
  345.    closenuff = delmin >> abs(periodicitycheck); /* for periodicity checking */
  346.    closenuff /= fudge;
  347.    rqlim_save = rqlim;
  348.    rqlim2 = sqrt(rqlim);
  349.    if (integerfractal)        /* for integer routines (lambda) */
  350.    {
  351.       lparm.x = parm.x * fudge;    /* real portion of Lambda */
  352.       lparm.y = parm.y * fudge;    /* imaginary portion of Lambda */
  353.       lparm2.x = parm2.x * fudge;  /* real portion of Lambda2 */
  354.       lparm2.y = parm2.y * fudge;  /* imaginary portion of Lambda2 */
  355.       llimit = rqlim * fudge;       /* stop if magnitude exceeds this */
  356.       if (llimit <= 0) llimit = 0x7fffffff; /* klooge for integer math */
  357.       llimit2 = rqlim2 * fudge;    /* stop if magnitude exceeds this */
  358.       lclosenuff = closenuff * fudge;    /* "close enough" value */
  359.       l16triglim = 8L<<16;       /* domain limit of fast trig functions */
  360.       linitorbit.x = initorbit.x * fudge;
  361.       linitorbit.y = initorbit.y * fudge;
  362.    }
  363.    resuming = (calc_status == 2);
  364.    if (!resuming) /* free resume_info memory if any is hanging around */
  365.    {
  366.       end_resume();
  367.       if (resave_flag) {
  368.      updatesavename(savename); /* do the pending increment */
  369.      resave_flag = started_resaves = 0;
  370.      }
  371.       calctime = 0;
  372.    }
  373.  
  374.    if (curfractalspecific->calctype != StandardFractal
  375.        && curfractalspecific->calctype != calcmand
  376.        && curfractalspecific->calctype != calcmandfp)
  377.    {
  378.       calctype = curfractalspecific->calctype; /* per_image can override */
  379.       symmetry = curfractalspecific->symmetry; /*   calctype & symmetry  */
  380.       plot = putcolor; /* defaults when setsymmetry not called or does nothing */
  381.       iystart = ixstart = yystart = xxstart = yybegin = 0;
  382.       iystop = yystop = ydots -1;
  383.       ixstop = xxstop = xdots -1;
  384.       calc_status = 1; /* mark as in-progress */
  385.       distest = 0; /* only standard escape time engine supports distest */
  386.       /* per_image routine is run here */
  387.       if (curfractalspecific->per_image())
  388.       { /* not a stand-alone */
  389.      /* next two lines in case periodicity changed */
  390.      closenuff = delmin >> abs(periodicitycheck); /* for periodicity checking */
  391.      closenuff /= fudge;
  392.      lclosenuff = closenuff * fudge;    /* "close enough" value */
  393.      setsymmetry(symmetry,0);
  394.      timer(0,calctype); /* non-standard fractal engine */
  395.       }
  396.       if (check_key())
  397.       {
  398.      if (calc_status == 1) /* calctype didn't set this itself, */
  399.         calc_status = 3;   /* so mark it interrupted, non-resumable */
  400.       }
  401.       else
  402.      calc_status = 4; /* no key, so assume it completed */
  403.    }
  404.    else /* standard escape-time engine */
  405.       timer(0,(int (*)())perform_worklist);
  406.    calctime += timer_interval;
  407.  
  408.    if(LogTable)
  409.    {
  410.       farmemfree(LogTable);
  411.       LogTable = NULL;
  412.    }
  413.    if(typespecific_workarea)
  414.    {
  415.       free_workarea();
  416.    }
  417.  
  418.    EXIT_OVLY;
  419.    return((calc_status == 4) ? 0 : -1);
  420. }
  421.  
  422. /**************** general escape-time engine routines *********************/
  423.  
  424. static void perform_worklist()
  425. {
  426.    int i;
  427.    long tmplong; /* this temp must be signed */
  428.  
  429.    if (potflag && pot16bit)
  430.    {
  431.       stdcalcmode = '1'; /* force 1 pass */
  432.       if (resuming == 0)
  433.      pot_startdisk();
  434.    }
  435.    if (stdcalcmode == 'b' && (curfractalspecific->flags & NOTRACE))
  436.       stdcalcmode = '1';
  437.    if (stdcalcmode == 'g' && (curfractalspecific->flags & NOGUESS))
  438.       stdcalcmode = '1';
  439.  
  440.    /* default setup a new worklist */
  441.    num_worklist = 1;
  442.    worklist[0].xxstart = 0;
  443.    worklist[0].yystart = worklist[0].yybegin = 0;
  444.    worklist[0].xxstop = xdots - 1;
  445.    worklist[0].yystop = ydots - 1;
  446.    worklist[0].pass = worklist[0].sym = 0;
  447.    if (resuming) /* restore worklist, if we can't the above will stay in place */
  448.    {
  449.       start_resume();
  450.       get_resume(sizeof(int),&num_worklist,sizeof(worklist),worklist,0);
  451.       end_resume();
  452.    }
  453.  
  454.    if (distest) /* setup stuff for distance estimator */
  455.    {
  456.       double ftemp,ftemp2;
  457.       dem_delta = sqr(delxx) + sqr(delyy2);
  458.       if ((ftemp = sqr(delyy) + sqr(delxx2)) > dem_delta)
  459.      dem_delta = ftemp;
  460.       if (distestwidth == 0)
  461.      distestwidth = 71;
  462.       ftemp = distestwidth;
  463.       dem_delta *= sqr(ftemp)/10000; /* multiply by thickness desired */
  464.       dem_width = ( sqrt( sqr(xxmax-xxmin) + sqr(xx3rd-xxmin) ) * ydots/xdots
  465.       + sqrt( sqr(yymax-yymin) + sqr(yy3rd-yymin) ) ) / distest;
  466.       ftemp = (rqlim < DEM_BAILOUT) ? DEM_BAILOUT : rqlim;
  467.       ftemp += 3; /* bailout plus just a bit */
  468.       ftemp2 = log(ftemp);
  469.       dem_toobig = sqr(ftemp) * sqr(ftemp2) * 4 / dem_delta;
  470.    }
  471.  
  472.    while (num_worklist > 0)
  473.    {
  474.       calctype = curfractalspecific->calctype; /* per_image can override */
  475.       symmetry = curfractalspecific->symmetry; /*   calctype & symmetry  */
  476.       plot = putcolor; /* defaults when setsymmetry not called or does nothing */
  477.  
  478.       /* pull top entry off worklist */
  479.       ixstart = xxstart = worklist[0].xxstart;
  480.       ixstop  = xxstop    = worklist[0].xxstop;
  481.       iystart = yystart = worklist[0].yystart;
  482.       iystop  = yystop    = worklist[0].yystop;
  483.       yybegin  = worklist[0].yybegin;
  484.       workpass = worklist[0].pass;
  485.       worksym  = worklist[0].sym;
  486.       --num_worklist;
  487.       for (i=0; i<num_worklist; ++i)
  488.      worklist[i] = worklist[i+1];
  489.  
  490.       calc_status = 1; /* mark as in-progress */
  491.  
  492.       curfractalspecific->per_image();
  493.  
  494.       /* some common initialization for escape-time pixel level routines */
  495.       closenuff = delmin >> abs(periodicitycheck); /* for periodicity checking */
  496.       closenuff /= fudge;
  497.       lclosenuff = closenuff * fudge;    /* "close enough" value */
  498.       kbdcount=max_kbdcount;
  499.       /* savedmask is for calcmand's periodicity checking */
  500.       savedmask = 0xC0000000; /* top 2 bits on */
  501.       tmplong = (delmin >> abs(periodicitycheck)) | 1;
  502.       while (tmplong > 0) /* while top bit not on */
  503.       {
  504.      tmplong <<= 1;
  505.      savedmask = (savedmask >> 1) | 0x80000000;
  506.       }
  507.  
  508.       setsymmetry(symmetry,1);
  509.  
  510.       /* call the appropriate escape-time engine */
  511.       switch (stdcalcmode)
  512.       {
  513.       case 't':
  514.      tesseral();
  515.      break;
  516.       case 'b':
  517.      bound_trace_main();
  518.      break;
  519.       case 'g':
  520.      solidguess();
  521.      break;
  522.       default:
  523.      OneOrTwoPass();
  524.       }
  525.  
  526.       if (check_key()) /* interrupted? */
  527.      break;
  528.    }
  529.  
  530.    if (num_worklist > 0)
  531.    {  /* interrupted, resumable */
  532.       alloc_resume(sizeof(worklist)+10,1);
  533.       put_resume(sizeof(int),&num_worklist,sizeof(worklist),worklist,0);
  534.    }
  535.    else
  536.       calc_status = 4; /* completed */
  537. }
  538.  
  539. static int OneOrTwoPass()
  540. {
  541.    int i;
  542.    totpasses = 1;
  543.    if (stdcalcmode == '2') totpasses = 2;
  544.    if (stdcalcmode == '2' && workpass == 0) /* do 1st pass of two */
  545.    {
  546.       if (StandardCalc(1) == -1)
  547.       {
  548.      add_worklist(xxstart,xxstop,yystart,yystop,row,0,worksym);
  549.      return(-1);
  550.       }
  551.       if (num_worklist > 0) /* worklist not empty, defer 2nd pass */
  552.       {
  553.      add_worklist(xxstart,xxstop,yystart,yystop,yystart,1,worksym);
  554.      return(0);
  555.       }
  556.       workpass = 1;
  557.       yybegin = yystart;
  558.    }
  559.    /* second or only pass */
  560.    if (StandardCalc(2) == -1)
  561.    {
  562.       i = yystop;
  563.       if (iystop != yystop) /* must be due to symmetry */
  564.      i -= row - iystart;
  565.       add_worklist(xxstart,xxstop,row,i,row,workpass,worksym);
  566.       return(-1);
  567.    }
  568.    return(0);
  569. }
  570.  
  571. static int _fastcall StandardCalc(int passnum)
  572. {
  573.    got_status = 0;
  574.    curpass = passnum;
  575.    row = yybegin;
  576.    while (row <= iystop)
  577.    {
  578.       currow = row;
  579.       reset_periodicity = 1;
  580.       col = ixstart;
  581.       while (col <= ixstop)
  582.       {
  583.          if(showdot>0)
  584.             (*plot) (col, row, showdot&(colors-1));
  585.  
  586.      /* on 2nd pass of two, skip even pts */
  587.      if (passnum == 1 || stdcalcmode == '1' || (row&1) != 0 || (col&1) != 0)
  588.      {
  589.         if ((*calctype)() == -1) /* StandardFractal(), calcmand() or calcmandfp() */
  590.            return(-1); /* interrupted */
  591.         reset_periodicity = 0;
  592.         if (passnum == 1) /* first pass, copy pixel and bump col */
  593.         {
  594.            if ((row&1) == 0 && row < iystop)
  595.            {
  596.           (*plot)(col,row+1,color);
  597.           if ((col&1) == 0 && col < ixstop)
  598.              (*plot)(col+1,row+1,color);
  599.            }
  600.            if ((col&1) == 0 && col < ixstop)
  601.           (*plot)(++col,row,color);
  602.         }
  603.      }
  604.      ++col;
  605.       }
  606.       if (passnum == 1 && (row&1) == 0)
  607.      ++row;
  608.       ++row;
  609.    }
  610.    return(0);
  611. }
  612.  
  613. int calcmand()        /* fast per pixel 1/2/b/g, called with row & col set */
  614. {
  615.    /* setup values from far array to avoid using es reg in calcmand.asm */
  616.    linitx = lx0[col] + lx1[row];
  617.    linity = ly0[row] + ly1[col];
  618.    if (calcmandasm() >= 0)
  619.    {
  620.       if (LogTable /* map color, but not if maxit & adjusted for inside,etc */
  621.       && (realcolor < maxit || (inside < 0 && color == maxit)))
  622.      color = LogTable[min(color, maxit)];
  623.       if (color >= colors) /* don't use color 0 unless from inside/outside */
  624.      if (colors < 16)
  625.         color &= andcolor;
  626.      else
  627.         color = ((color - 1) % andcolor) + 1;  /* skip color zero */
  628.       (*plot) (col, row, color);
  629.    }
  630.    return (color);
  631. }
  632.  
  633. /************************************************************************/
  634. /* added by Wes Loewer - sort of a floating point version of calcmand() */
  635. /* can also handle invert, any rqlim, potflag, zmag, epsilon cross,     */
  636. /* and all the current outside options    -Wes Loewer 11/03/91          */
  637. /************************************************************************/
  638. int calcmandfp()
  639. {
  640.    if(invert)
  641.       invertz2(&init);
  642.    else
  643.    {
  644.       init.y = dy0[row]+dy1[col];
  645.       init.x = dx0[col]+dx1[row];
  646.    }
  647.    if (calcmandfpasm() >= 0)
  648.    {
  649.       if (potflag)
  650.      color = potential(magnitude, realcolor);
  651.       if (LogTable /* map color, but not if maxit & adjusted for inside,etc */
  652.           && (realcolor < maxit || (inside < 0 && color == maxit)))
  653.         color = LogTable[min(color, maxit)];
  654.      if (color >= colors) /* don't use color 0 unless from inside/outside */
  655.         if (colors < 16)
  656.            color &= andcolor;
  657.         else
  658.            color = ((color - 1) % andcolor) + 1;  /* skip color zero */
  659.       (*plot) (col, row, color);
  660.    }
  661.    return (color);
  662. }
  663.  
  664. #define green 2
  665. #define yellow 6
  666. int StandardFractal()    /* per pixel 1/2/b/g, called with row & col set */
  667. {
  668.    double tantable[16];
  669.    int hooper;
  670.    double close;
  671.    long lclose;
  672.    int cyclelen = -1;
  673.    int savedcolor;
  674.    int caught_a_cycle;
  675.    int savedand, savedincr;    /* for periodicity checking */
  676.    struct lcomplex lsaved;
  677.    int i, attracted;
  678.    struct lcomplex lat;
  679.    struct complex  at;
  680.    struct complex deriv;
  681.    int dem_color;
  682.    struct complex dem_new;
  683.    close = .01;
  684.    lclose = close*fudge;
  685.    if(inside == -101)
  686.    {
  687.       int i;
  688.       for(i=0;i<16;i++)
  689.      tantable[i] = 0.0;
  690.    }
  691.    else if(inside == -100)
  692.    {
  693.       close = .01;
  694.       lclose = close*fudge;
  695.    }
  696.    if (periodicitycheck == 0 || inside == -59 || inside == -101)
  697.       oldcolor = 32767;     /* don't check periodicity at all */
  698.    else if (inside == -102)    /* for display-periodicity */
  699.       oldcolor = maxit*4/5;    /* don't check until nearly done */
  700.    else if (reset_periodicity)
  701.       oldcolor = 250;        /* don't check periodicity 1st 250 iterations */
  702.  
  703.    /* really fractal specific, but we'll leave it here */
  704.    if (!integerfractal)
  705.    {
  706.       if (useinitorbit == 1)
  707.      saved = initorbit;
  708.       else {
  709.      saved.x = 0;
  710.      saved.y = 0;
  711.       }
  712.       init.y = dy0[row] + dy1[col];
  713.       if (distest)
  714.       {
  715.      rqlim = rqlim_save;          /* start with regular bailout */
  716.      if (distest != 1 || colors == 2) /* not doing regular outside colors */
  717.         if (rqlim < DEM_BAILOUT)      /* so go straight for dem bailout */
  718.            rqlim = DEM_BAILOUT;
  719.      deriv.x = 1;
  720.      deriv.y = 0;
  721.      magnitude = 0;
  722.      dem_color = -1;
  723.       }
  724.    }
  725.    else
  726.    {
  727.       if (useinitorbit == 1)
  728.      lsaved = linitorbit;
  729.       else {
  730.      lsaved.x = 0;
  731.      lsaved.y = 0;
  732.       }
  733.       linit.y = ly0[row] + ly1[col];
  734.    }
  735.    orbit_ptr = 0;
  736.    color = 0;
  737.    if(fractype==JULIAFP || fractype==JULIA)
  738.       color = -1;
  739.    caught_a_cycle = 0;
  740.    if (inside == -102) {
  741.        savedand = 16;        /* begin checking every 16th cycle */
  742.    } else {
  743.        savedand = 1;        /* begin checking every other cycle */
  744.    }
  745.    savedincr = 1;        /* start checking the very first time */
  746.  
  747.    if (inside <= -60 && inside >= -61)
  748.    {
  749.       magnitude = lmagnitud = 0;
  750.       min_orbit = 100000.0;
  751.    }
  752.    overflow = 0;        /* reset integer math overflow flag */
  753.  
  754.    curfractalspecific->per_pixel(); /* initialize the calculations */
  755.  
  756.    attracted = FALSE;
  757.    while (++color < maxit)
  758.    {
  759.       if(showdot>0)
  760.          (*plot) (col, row, showdot&(colors-1));
  761.  
  762.       /* calculation of one orbit goes here */
  763.       /* input in "old" -- output in "new" */
  764.  
  765.       if (distest)
  766.       {
  767.      double ftemp;
  768.      /* Distance estimator for points near Mandelbrot set */
  769.      /* Original code by Phil Wilson, hacked around by PB */
  770.      /* Algorithms from Peitgen & Saupe, Science of Fractal Images, p.198 */
  771.      ftemp     = 2 * (old.x * deriv.x - old.y * deriv.y) + 1;
  772.      deriv.y = 2 * (old.y * deriv.x + old.x * deriv.y);
  773.      deriv.x = ftemp;
  774.      if (sqr(deriv.x)+sqr(deriv.y) > dem_toobig)
  775.         break;
  776.      /* if above exit taken, the later test vs dem_delta will place this
  777.             point on the boundary, because mag(old)<bailout just now */
  778.      if (curfractalspecific->orbitcalc())
  779.      {
  780.         if (dem_color < 0)          /* note "regular" color for later */
  781.         {
  782.            dem_color = color;
  783.            dem_new = new;
  784.         }
  785.         if (rqlim >= DEM_BAILOUT  /* exit if past real bailout */
  786.         || magnitude >= (rqlim = DEM_BAILOUT) /* reset to real bailout */
  787.         || magnitude == 0)     /* exit if type doesn't "floatbailout" */
  788.            break;
  789.         old = new;              /* carry on till past real bailout */
  790.      }
  791.       }
  792.       else /* the usual case */
  793.      if (curfractalspecific->orbitcalc() && inside != -101)
  794.         break;
  795.       if (show_orbit)
  796.      if (!integerfractal)
  797.         plot_orbit(new.x, new.y, -1);
  798.      else
  799.         iplot_orbit(lnew.x, lnew.y, -1);
  800.       if(inside == -101)
  801.       {
  802.      if(0 < color && color < 16)
  803.      {
  804.      if (integerfractal)
  805.      {
  806.            new.x = lnew.x;
  807.            new.x /= fudge;
  808.            new.y = lnew.y;
  809.            new.y /= fudge;
  810.         }
  811.         tantable[color-1] = new.y/(new.x+.000001);
  812.      }
  813.       }
  814.       else if(inside == -100)
  815.       {
  816.      hooper = 0;
  817.      if(integerfractal)
  818.      {
  819.         if(labs(lnew.x) < lclose)
  820.         {
  821.            hooper = 1; /* close to y axis */
  822.            goto plot_inside;
  823.         }
  824.         else if(labs(lnew.y) < lclose)
  825.         {
  826.            hooper = 2; /* close to x axis */
  827.            goto plot_inside;
  828.         }
  829.      }
  830.      else
  831.      {
  832.         if(fabs(new.x) < close)
  833.         {
  834.            hooper = 1; /* close to y axis */
  835.            goto plot_inside;
  836.         }
  837.         else if(fabs(new.y) < close)
  838.         {
  839.            hooper = 2; /* close to x axis */
  840.            goto plot_inside;
  841.         }
  842.      }
  843.       }
  844.       else if (inside <= -60 && inside >= -61)
  845.       {
  846.      if (integerfractal)
  847.      {
  848.         if (lmagnitud == 0)
  849.            lmagnitud = lsqr(lnew.x) + lsqr(lnew.y);
  850.         magnitude = lmagnitud;
  851.         magnitude = magnitude / fudge;
  852.      }
  853.      else
  854.         if (magnitude == 0.0)
  855.            magnitude = sqr(new.x) + sqr(new.y);
  856.      if (magnitude < min_orbit)
  857.      {
  858.         min_orbit = magnitude;
  859.         min_index = color + 1;
  860.      }
  861.       }
  862.  
  863.  
  864.       if (attractors > 0)    /* finite attractor in the list   */
  865.       {             /* NOTE: Integer code is UNTESTED */
  866.      if (integerfractal)
  867.      {
  868.         for (i = 0; i < attractors; i++)
  869.         {
  870.         lat.x = lnew.x - lattr[i].x;
  871.         lat.x = lsqr(lat.x);
  872.         if (lat.x < l_at_rad)
  873.         {
  874.            lat.y = lnew.y - lattr[i].y;
  875.            lat.y = lsqr(lat.y);
  876.            if (lat.y < l_at_rad)
  877.            {
  878.               if ((lat.x + lat.y) < l_at_rad)
  879.               {
  880.              attracted = TRUE;
  881.              if (finattract<0) color = (color%attrperiod[i])+1;
  882.              break;
  883.               }
  884.            }
  885.         }
  886.         }
  887.      }
  888.      else
  889.      {
  890.         for (i = 0; i < attractors; i++)
  891.         {
  892.         at.x = new.x - attr[i].x;
  893.         at.x = sqr(at.x);
  894.         if (at.x < f_at_rad)
  895.         {
  896.            at.y = new.y - attr[i].y;
  897.            at.y = sqr(at.y);
  898.            if ( at.y < f_at_rad)
  899.            {
  900.               if ((at.x + at.y) < f_at_rad)
  901.               {
  902.              attracted = TRUE;
  903.              if (finattract<0) color = (color%attrperiod[i])+1;
  904.              break;
  905.               }
  906.            }
  907.         }
  908.         }
  909.      }
  910.      if (attracted)
  911.         break;        /* AHA! Eaten by an attractor */
  912.       }
  913.  
  914.       if (color > oldcolor)    /* check periodicity */
  915.       {
  916.      if ((color & savedand) == 0)         /* time to save a new value */
  917.      {
  918.         savedcolor = color;
  919.         if (!integerfractal)
  920.            saved = new;  /* floating pt fractals */
  921.         else
  922.            lsaved = lnew;/* integer fractals */
  923.         if (--savedincr == 0)    /* time to lengthen the periodicity? */
  924.         {
  925.            savedand = (savedand << 1) + 1;         /* longer periodicity */
  926.            savedincr = 4;/* restart counter */
  927.         }
  928.      }
  929.      else             /* check against an old save */
  930.      {
  931.         if (!integerfractal)     /* floating-pt periodicity chk */
  932.         {
  933.            if (fabs(saved.x - new.x) < closenuff)
  934.           if (fabs(saved.y - new.y) < closenuff)
  935.           {
  936.              caught_a_cycle = 1;
  937.              cyclelen = color-savedcolor;
  938.              color = maxit - 1;
  939.           }
  940.         }
  941.         else         /* integer periodicity check */
  942.         {
  943.            if (labs(lsaved.x - lnew.x) < lclosenuff)
  944.           if (labs(lsaved.y - lnew.y) < lclosenuff)
  945.           {
  946.              caught_a_cycle = 1;
  947.              cyclelen = color-savedcolor;
  948.              color = maxit - 1;
  949.           }
  950.         }
  951.      }
  952.       }
  953.    }
  954.  
  955.    if (show_orbit)
  956.       scrub_orbit();
  957.  
  958.    realcolor = color;        /* save this before we start adjusting it */
  959.    if (color >= maxit)
  960.       oldcolor = 0;        /* check periodicity immediately next time */
  961.    else
  962.    {
  963.       oldcolor = color + 10;    /* check when past this + 10 next time */
  964.       if (color == 0)
  965.      color = 1;        /* needed to make same as calcmand */
  966.    }
  967.  
  968.    if (potflag)
  969.    {
  970.       if (integerfractal)    /* adjust integer fractals */
  971.       {
  972.      new.x = ((double)lnew.x) / fudge;
  973.      new.y = ((double)lnew.y) / fudge;
  974.       }
  975.       magnitude = sqr(new.x) + sqr(new.y);
  976.       color = potential(magnitude, color);
  977.       if (LogTable)
  978.      color = LogTable[min(color, maxit)];
  979.       goto plot_pixel;        /* skip any other adjustments */
  980.    }
  981.  
  982.    if (color >= maxit)        /* an "inside" point */
  983.       goto plot_inside;     /* distest, decomp, biomorph don't apply */
  984.  
  985.  
  986.    if (outside < -1)  /* these options by Richard Hughes modified by TW */
  987.    {
  988.       if (integerfractal)
  989.       {
  990.      new.x = ((double)lnew.x) / fudge;
  991.      new.y = ((double)lnew.y) / fudge;
  992.       }
  993.       /* Add 7 to overcome negative values on the MANDEL    */
  994.       if (outside == -2)         /* "real" */
  995.      color += new.x + 7;
  996.       else if (outside == -3)         /* "imag" */
  997.      color += new.y + 7;
  998.       else if (outside == -4  && new.y)  /* "mult" */
  999.       color *= (new.x/new.y);
  1000.       else if (outside == -5)         /* "sum" */
  1001.       color += new.x + new.y;
  1002.  
  1003.       /* eliminate negative colors & wrap arounds */
  1004.       if (color < 0 || color > maxit)
  1005.       color = 0;
  1006.    }
  1007.  
  1008.    if (distest)
  1009.    {
  1010.       double dist,temp;
  1011.       dist = sqr(new.x) + sqr(new.y);
  1012.       temp = log(dist);
  1013.       dist = dist * sqr(temp) / ( sqr(deriv.x) + sqr(deriv.y) );
  1014.       if (dist < dem_delta)    /* point is on the edge */
  1015.       {
  1016.      if (distest > 0)
  1017.         goto plot_inside;    /* show it as an inside point */
  1018.      color = 0 - distest;    /* show boundary as specified color */
  1019.      goto plot_pixel;    /* no further adjustments apply */
  1020.       }
  1021.       if (colors == 2)
  1022.       {
  1023.      color = !inside;    /* the only useful distest 2 color use */
  1024.      goto plot_pixel;    /* no further adjustments apply */
  1025.       }
  1026.       if (distest > 1)        /* pick color based on distance */
  1027.       {
  1028.      color = sqrt(dist / dem_width + 1);
  1029.      goto plot_pixel;    /* no further adjustments apply */
  1030.       }
  1031.       color = dem_color;    /* use pixel's "regular" color */
  1032.       new = dem_new;
  1033.    }
  1034.  
  1035.    if (decomp[0] > 0)
  1036.       decomposition();
  1037.    else if (biomorph != -1)
  1038.    {
  1039.       if (integerfractal)
  1040.       {
  1041.      if (labs(lnew.x) < llimit2 || labs(lnew.y) < llimit2)
  1042.         color = biomorph;
  1043.       }
  1044.       else
  1045.      if (fabs(new.x) < rqlim2 || fabs(new.y) < rqlim2)
  1046.         color = biomorph;
  1047.    }
  1048.  
  1049.    if (outside >= 0 && attracted == FALSE) /* merge escape-time stripes */
  1050.       color = outside;
  1051.    else if (LogTable)
  1052.       color = LogTable[min(color, maxit)];
  1053.    goto plot_pixel;
  1054.  
  1055.    plot_inside: /* we're "inside" */
  1056.    if (periodicitycheck < 0 && caught_a_cycle)
  1057.       color = 7;           /* show periodicity */
  1058.    else if (inside >= 0)
  1059.       color = inside;           /* set to specified color, ignore logpal */
  1060.    else
  1061.    {
  1062.       if(inside == -101)
  1063.       {
  1064.      int i;
  1065.      double diff;
  1066.      color = 0;
  1067.      for(i=1;i<16;i++)
  1068.      {
  1069.         diff = tantable[0] - tantable[i];
  1070.         if(fabs(diff) < .05)
  1071.         {
  1072.            color = i;
  1073.            break;
  1074.         }
  1075.      }
  1076.       }
  1077.       else if(inside== -102) {
  1078.       if (cyclelen>0) {
  1079.           color = cyclelen;
  1080.       } else {
  1081.           color = maxit;
  1082.       }
  1083.       }
  1084.       else if(inside == -100)
  1085.       {
  1086.      if(hooper==1)
  1087.         color = green;
  1088.      else if(hooper==2)
  1089.         color = yellow;
  1090.      else
  1091.         color = maxit;
  1092.      if (show_orbit)
  1093.         scrub_orbit();
  1094.       }
  1095.       else if (inside == -60)
  1096.      color = sqrt(min_orbit) * 75;
  1097.       else if (inside == -61)
  1098.      color = min_index;
  1099.       else if (inside == -59)
  1100.       {
  1101.      if (integerfractal)
  1102.      {
  1103.         /*
  1104.         new.x = ((double)lnew.x) / fudge;
  1105.         new.y = ((double)lnew.y) / fudge;
  1106.         */
  1107.         color = (((double)lsqr(lnew.x))/fudge + ((double)lsqr(lnew.y))/fudge) * (maxit>>1) + 1;
  1108.      }
  1109.      else
  1110.         color = (sqr(new.x) + sqr(new.y)) * (maxit>>1) + 1;
  1111.       }
  1112.       else /* inside == -1 */
  1113.      color = maxit;
  1114.       if (LogTable)
  1115.      color = LogTable[min(color, maxit)];
  1116.    }
  1117.  
  1118.    plot_pixel:
  1119.  
  1120.    if (color >= colors) /* don't use color 0 unless from inside/outside */
  1121.       if (colors < 16)
  1122.      color &= andcolor;
  1123.       else
  1124.      color = ((color - 1) % andcolor) + 1;    /* skip color zero */
  1125.  
  1126.    (*plot) (col, row, color);
  1127.  
  1128.    if ((kbdcount -= realcolor) <= 0)
  1129.    {
  1130.       if (check_key())
  1131.      return (-1);
  1132.       kbdcount = max_kbdcount;
  1133.    }
  1134.  
  1135.    return (color);
  1136. }
  1137. #undef green
  1138. #undef yellow
  1139.  
  1140.  
  1141.  
  1142. /**************** standardfractal doodad subroutines *********************/
  1143.  
  1144. static void decomposition()
  1145. {
  1146.    static double cos45       = 0.70710678118654750; /* cos 45    degrees */
  1147.    static double sin45       = 0.70710678118654750; /* sin 45    degrees */
  1148.    static double cos22_5   = 0.92387953251128670; /* cos 22.5    degrees */
  1149.    static double sin22_5   = 0.38268343236508980; /* sin 22.5    degrees */
  1150.    static double cos11_25  = 0.98078528040323040; /* cos 11.25    degrees */
  1151.    static double sin11_25  = 0.19509032201612820; /* sin 11.25    degrees */
  1152.    static double cos5_625  = 0.99518472667219690; /* cos 5.625    degrees */
  1153.    static double sin5_625  = 0.09801714032956060; /* sin 5.625    degrees */
  1154.    static double tan22_5   = 0.41421356237309500; /* tan 22.5    degrees */
  1155.    static double tan11_25  = 0.19891236737965800; /* tan 11.25    degrees */
  1156.    static double tan5_625  = 0.09849140335716425; /* tan 5.625    degrees */
  1157.    static double tan2_8125 = 0.04912684976946725; /* tan 2.8125 degrees */
  1158.    static double tan1_4063 = 0.02454862210892544; /* tan 1.4063 degrees */
  1159.    static long lcos45      ; /* cos 45      degrees */
  1160.    static long lsin45      ; /* sin 45      degrees */
  1161.    static long lcos22_5   ; /* cos 22.5   degrees */
  1162.    static long lsin22_5   ; /* sin 22.5   degrees */
  1163.    static long lcos11_25  ; /* cos 11.25  degrees */
  1164.    static long lsin11_25  ; /* sin 11.25  degrees */
  1165.    static long lcos5_625  ; /* cos 5.625  degrees */
  1166.    static long lsin5_625  ; /* sin 5.625  degrees */
  1167.    static long ltan22_5   ; /* tan 22.5   degrees */
  1168.    static long ltan11_25  ; /* tan 11.25  degrees */
  1169.    static long ltan5_625  ; /* tan 5.625  degrees */
  1170.    static long ltan2_8125 ; /* tan 2.8125 degrees */
  1171.    static long ltan1_4063 ; /* tan 1.4063 degrees */
  1172.    static reset_fudge = -1;
  1173.    int temp = 0;
  1174.    int i;
  1175.    struct lcomplex lalt;
  1176.    struct complex alt;
  1177.    color = 0;
  1178.    if (integerfractal) /* the only case */
  1179.    {
  1180.       if (reset_fudge != fudge)
  1181.       {
  1182.      reset_fudge = fudge;
  1183.      lcos45     = cos45     *fudge;
  1184.      lsin45     = sin45     *fudge;
  1185.      lcos22_5   = cos22_5     *fudge;
  1186.      lsin22_5   = sin22_5     *fudge;
  1187.      lcos11_25  = cos11_25     *fudge;
  1188.      lsin11_25  = sin11_25     *fudge;
  1189.      lcos5_625  = cos5_625     *fudge;
  1190.      lsin5_625  = sin5_625     *fudge;
  1191.      ltan22_5   = tan22_5     *fudge;
  1192.      ltan11_25  = tan11_25     *fudge;
  1193.      ltan5_625  = tan5_625     *fudge;
  1194.      ltan2_8125 = tan2_8125  *fudge;
  1195.      ltan1_4063 = tan1_4063  *fudge;
  1196.       }
  1197.       if (lnew.y < 0)
  1198.       {
  1199.      temp = 2;
  1200.      lnew.y = -lnew.y;
  1201.       }
  1202.  
  1203.       if (lnew.x < 0)
  1204.       {
  1205.      ++temp;
  1206.      lnew.x = -lnew.x;
  1207.       }
  1208.  
  1209.       if (decomp[0] >= 8)
  1210.       {
  1211.      temp <<= 1;
  1212.      if (lnew.x < lnew.y)
  1213.      {
  1214.         ++temp;
  1215.         lalt.x = lnew.x; /* just */
  1216.         lnew.x = lnew.y; /* swap */
  1217.         lnew.y = lalt.x; /* them */
  1218.      }
  1219.  
  1220.      if (decomp[0] >= 16)
  1221.      {
  1222.         temp <<= 1;
  1223.         if (multiply(lnew.x,ltan22_5,bitshift) < lnew.y)
  1224.         {
  1225.            ++temp;
  1226.            lalt = lnew;
  1227.            lnew.x = multiply(lalt.x,lcos45,bitshift) +
  1228.            multiply(lalt.y,lsin45,bitshift);
  1229.            lnew.y = multiply(lalt.x,lsin45,bitshift) -
  1230.            multiply(lalt.y,lcos45,bitshift);
  1231.         }
  1232.  
  1233.         if (decomp[0] >= 32)
  1234.         {
  1235.            temp <<= 1;
  1236.            if (multiply(lnew.x,ltan11_25,bitshift) < lnew.y)
  1237.            {
  1238.           ++temp;
  1239.           lalt = lnew;
  1240.           lnew.x = multiply(lalt.x,lcos22_5,bitshift) +
  1241.               multiply(lalt.y,lsin22_5,bitshift);
  1242.           lnew.y = multiply(lalt.x,lsin22_5,bitshift) -
  1243.               multiply(lalt.y,lcos22_5,bitshift);
  1244.            }
  1245.  
  1246.            if (decomp[0] >= 64)
  1247.            {
  1248.           temp <<= 1;
  1249.           if (multiply(lnew.x,ltan5_625,bitshift) < lnew.y)
  1250.           {
  1251.              ++temp;
  1252.              lalt = lnew;
  1253.              lnew.x = multiply(lalt.x,lcos11_25,bitshift) +
  1254.              multiply(lalt.y,lsin11_25,bitshift);
  1255.              lnew.y = multiply(lalt.x,lsin11_25,bitshift) -
  1256.              multiply(lalt.y,lcos11_25,bitshift);
  1257.           }
  1258.  
  1259.           if (decomp[0] >= 128)
  1260.           {
  1261.              temp <<= 1;
  1262.              if (multiply(lnew.x,ltan2_8125,bitshift) < lnew.y)
  1263.              {
  1264.             ++temp;
  1265.             lalt = lnew;
  1266.             lnew.x = multiply(lalt.x,lcos5_625,bitshift) +
  1267.                 multiply(lalt.y,lsin5_625,bitshift);
  1268.             lnew.y = multiply(lalt.x,lsin5_625,bitshift) -
  1269.                 multiply(lalt.y,lcos5_625,bitshift);
  1270.              }
  1271.  
  1272.              if (decomp[0] == 256)
  1273.              {
  1274.             temp <<= 1;
  1275.             if (multiply(lnew.x,ltan1_4063,bitshift) < lnew.y)
  1276.                if ((lnew.x*ltan1_4063 < lnew.y))
  1277.                   ++temp;
  1278.              }
  1279.           }
  1280.            }
  1281.         }
  1282.      }
  1283.       }
  1284.    }
  1285.    else /* double case */
  1286.    {
  1287.       if (new.y < 0)
  1288.       {
  1289.      temp = 2;
  1290.      new.y = -new.y;
  1291.       }
  1292.       if (new.x < 0)
  1293.       {
  1294.      ++temp;
  1295.      new.x = -new.x;
  1296.       }
  1297.       if (decomp[0] >= 8)
  1298.       {
  1299.      temp <<= 1;
  1300.      if (new.x < new.y)
  1301.      {
  1302.         ++temp;
  1303.         alt.x = new.x; /* just */
  1304.         new.x = new.y; /* swap */
  1305.         new.y = alt.x; /* them */
  1306.      }
  1307.      if (decomp[0] >= 16)
  1308.      {
  1309.         temp <<= 1;
  1310.         if (new.x*tan22_5 < new.y)
  1311.         {
  1312.            ++temp;
  1313.            alt = new;
  1314.            new.x = alt.x*cos45 + alt.y*sin45;
  1315.            new.y = alt.x*sin45 - alt.y*cos45;
  1316.         }
  1317.  
  1318.         if (decomp[0] >= 32)
  1319.         {
  1320.            temp <<= 1;
  1321.            if (new.x*tan11_25 < new.y)
  1322.            {
  1323.           ++temp;
  1324.           alt = new;
  1325.           new.x = alt.x*cos22_5 + alt.y*sin22_5;
  1326.           new.y = alt.x*sin22_5 - alt.y*cos22_5;
  1327.            }
  1328.  
  1329.            if (decomp[0] >= 64)
  1330.            {
  1331.           temp <<= 1;
  1332.           if (new.x*tan5_625 < new.y)
  1333.           {
  1334.              ++temp;
  1335.              alt = new;
  1336.              new.x = alt.x*cos11_25 + alt.y*sin11_25;
  1337.              new.y = alt.x*sin11_25 - alt.y*cos11_25;
  1338.           }
  1339.  
  1340.           if (decomp[0] >= 128)
  1341.           {
  1342.              temp <<= 1;
  1343.              if (new.x*tan2_8125 < new.y)
  1344.              {
  1345.             ++temp;
  1346.             alt = new;
  1347.             new.x = alt.x*cos5_625 + alt.y*sin5_625;
  1348.             new.y = alt.x*sin5_625 - alt.y*cos5_625;
  1349.              }
  1350.  
  1351.              if (decomp[0] == 256)
  1352.              {
  1353.             temp <<= 1;
  1354.             if ((new.x*tan1_4063 < new.y))
  1355.                ++temp;
  1356.              }
  1357.           }
  1358.            }
  1359.         }
  1360.      }
  1361.       }
  1362.    }
  1363.    for (i = 1; temp > 0; ++i)
  1364.    {
  1365.       if (temp & 1)
  1366.      color = (1 << i) - 1 - color;
  1367.       temp >>= 1;
  1368.    }
  1369.    if (decomp[0] == 2)
  1370.       color &= 1;
  1371.    if (colors > decomp[0])
  1372.       color++;
  1373. }
  1374.  
  1375. /******************************************************************/
  1376. /* Continuous potential calculation for Mandelbrot and Julia      */
  1377. /* Reference: Science of Fractal Images p. 190.           */
  1378. /* Special thanks to Mark Peterson for his "MtMand" program that  */
  1379. /* beautifully approximates plate 25 (same reference) and spurred */
  1380. /* on the inclusion of similar capabilities in FRACTINT.      */
  1381. /*                                  */
  1382. /* The purpose of this function is to calculate a color value      */
  1383. /* for a fractal that varies continuously with the screen pixels  */
  1384. /* locations for better rendering in 3D.              */
  1385. /*                                  */
  1386. /* Here "magnitude" is the modulus of the orbit value at          */
  1387. /* "iterations". The potparms[] are user-entered paramters        */
  1388. /* controlling the level and slope of the continuous potential      */
  1389. /* surface. Returns color.  - Tim Wegner 6/25/89          */
  1390. /*                                  */
  1391. /*               -- Change history --              */
  1392. /*                                  */
  1393. /* 09/12/89   - added floatflag support and fixed float underflow */
  1394. /*                                  */
  1395. /******************************************************************/
  1396.  
  1397. static int _fastcall potential(double mag, int iterations)
  1398. {
  1399.    float f_mag,f_tmp,pot;
  1400.    double d_tmp;
  1401.    int i_pot;
  1402.    long l_pot;
  1403.    extern char floatflag;
  1404.  
  1405.    if(iterations < maxit)
  1406.    {
  1407.       pot = i_pot = iterations+2;
  1408.       if(i_pot <= 0 || mag <= 1.0)
  1409.      pot = 0.0;
  1410.       else /* pot = log(mag) / pow(2.0, (double)pot); */
  1411.       {
  1412.      if(i_pot < 120 && !floatflag) /* empirically determined limit of fShift */
  1413.      {
  1414.         f_mag = mag;
  1415.         fLog14(f_mag,f_tmp); /* this SHOULD be non-negative */
  1416.         fShift(f_tmp,(char)-i_pot,pot);
  1417.      }
  1418.      else
  1419.      {
  1420.         d_tmp = log(mag)/(double)pow(2.0,(double)pot);
  1421.         if(d_tmp > FLT_MIN) /* prevent float type underflow */
  1422.            pot = d_tmp;
  1423.         else
  1424.            pot = 0.0;
  1425.      }
  1426.       }
  1427.       /* following transformation strictly for aesthetic reasons */
  1428.       /* meaning of parameters:
  1429.         potparam[0] -- zero potential level - highest color -
  1430.         potparam[1] -- slope multiplier -- higher is steeper
  1431.         potparam[2] -- rqlim value if changeable (bailout for modulus) */
  1432.  
  1433.       if(pot > 0.0)
  1434.       {
  1435.      if(floatflag)
  1436.         pot = (float)sqrt((double)pot);
  1437.      else
  1438.      {
  1439.         fSqrt14(pot,f_tmp);
  1440.         pot = f_tmp;
  1441.      }
  1442.      pot = potparam[0] - pot*potparam[1] - 1.0;
  1443.       }
  1444.       else
  1445.      pot = potparam[0] - 1.0;
  1446.       if(pot < 1.0)
  1447.      pot = 1.0; /* avoid color 0 */
  1448.    }
  1449.    else if(inside >= 0)
  1450.       pot = inside;
  1451.    else /* inside < 0 implies inside=maxit, so use 1st pot param instead */
  1452.       pot = potparam[0];
  1453.  
  1454.    i_pot = (l_pot = pot * 256) >> 8;
  1455.    if(i_pot >= colors)
  1456.    {
  1457.       i_pot = colors - 1;
  1458.       l_pot = 255;
  1459.    }
  1460.  
  1461.    if(pot16bit)
  1462.    {
  1463.       if (dotmode != 11) /* if putcolor won't be doing it for us */
  1464.      writedisk(col+sxoffs,row+syoffs,i_pot);
  1465.       writedisk(col+sxoffs,row+sydots+syoffs,(int)l_pot);
  1466.    }
  1467.  
  1468.    return(i_pot);
  1469. }
  1470.  
  1471.  
  1472. /******************* boundary trace method ***************************
  1473. Fractint's original btm was written by David Guenther.  There were a few
  1474. rare circumstances in which the original btm would not trace or fill
  1475. correctly, even on Mandelbrot Sets.  The code below was adapted from
  1476. "Mandelbrot Sets by Wesley Loewer" (see calmanfp.asm) which was written
  1477. before I was introduced to Fractint.  It should be noted that without
  1478. David Guenther's implimentation of a btm, I doubt that I would have been
  1479. able to impliment my own code into Fractint.  There are several things in
  1480. the following code that are not original with me but came from David
  1481. Guenther's code.  I've noted these places with the initials DG.
  1482.  
  1483.                     Wesley Loewer 3/8/92
  1484. *********************************************************************/
  1485. #define bkcolor 0  /* I have some ideas for the future with this. -Wes */
  1486. #define advance_match()     coming_from = ((going_to = (going_to - 1) & 0x03) - 1) & 0x03
  1487. #define advance_no_match()  going_to = (going_to + 1) & 0x03
  1488.  
  1489. static
  1490. int  bound_trace_main()
  1491.     {
  1492.     enum direction coming_from;
  1493.     int match_found, continue_loop;
  1494.     int trail_color, fillcolor_used, last_fillcolor_used;
  1495.     int max_putline_length;
  1496.     int right, left, length;
  1497.  
  1498.     if (inside == 0 || outside == 0)
  1499.     {
  1500.     static char far msg[]=
  1501.         "Boundary tracing cannot be used with inside=0 or outside=0.";
  1502.     stopmsg(0,msg);
  1503.     return(-1);
  1504.     }
  1505.     if (colors < 16)
  1506.     {
  1507.     static char far msg[]=
  1508.         "Boundary tracing cannot be used with < 16 colors.";
  1509.     stopmsg(0,msg);
  1510.     return(-1);
  1511.     }
  1512.  
  1513.     got_status = 2;
  1514.     max_putline_length = 0; /* reset max_putline_length */
  1515.     for (currow = iystart; currow <= iystop; currow++)
  1516.     {
  1517.     reset_periodicity = 1; /* reset for a new row */
  1518.     color = bkcolor;
  1519.     for (curcol = ixstart; curcol <= ixstop; curcol++)
  1520.         {
  1521.         if (getcolor(curcol, currow) != bkcolor)
  1522.         continue;
  1523.  
  1524.         trail_color = color;
  1525.         row = currow;
  1526.         col = curcol;
  1527.         if ((*calctype)()== -1) /* color, row, col are global */
  1528.         {
  1529.         if (iystop != yystop)  /* DG */
  1530.            iystop = yystop - (currow - yystart); /* allow for sym */
  1531.         add_worklist(xxstart,xxstop,currow,iystop,currow,0,worksym);
  1532.         return -1;
  1533.         }
  1534.         reset_periodicity = 0; /* normal periodicity checking */
  1535.  
  1536.         /*
  1537.         This next line may cause a few more pixels to be calculated,
  1538.         but at the savings of quite a bit of overhead
  1539.         */
  1540.         if (color != trail_color)  /* DG */
  1541.         continue;
  1542.  
  1543.         /* sweep clockwise to trace outline */
  1544.         trail_row = currow;
  1545.         trail_col = curcol;
  1546.         trail_color = color;
  1547.         fillcolor_used = fillcolor > 0 ? fillcolor : trail_color;
  1548.         coming_from = West;
  1549.         going_to = East;
  1550.         match_found = 0;
  1551.         continue_loop = TRUE;
  1552.         do
  1553.         {
  1554.         step_col_row();
  1555.         if (row >= currow
  1556.             && col >= ixstart
  1557.             && col <= ixstop
  1558.             && row <= iystop)
  1559.             {
  1560.             /* the order of operations in this next line is critical */
  1561.             if ((color = getcolor(col, row)) == bkcolor && (*calctype)()== -1)
  1562.                 /* color, row, col are global for (*calctype)() */
  1563.             {
  1564.             if (iystop != yystop)  /* DG */
  1565.                iystop = yystop - (currow - yystart); /* allow for sym */
  1566.             add_worklist(xxstart,xxstop,currow,iystop,currow,0,worksym);
  1567.             return -1;
  1568.             }
  1569.             else if (color == trail_color)
  1570.             {
  1571.             match_found++;
  1572.             trail_row = row;
  1573.             trail_col = col;
  1574.             advance_match();
  1575.             }
  1576.             else
  1577.             {
  1578.             advance_no_match();
  1579.             continue_loop = going_to != coming_from || match_found;
  1580.             }
  1581.             }
  1582.         else
  1583.             {
  1584.             advance_no_match();
  1585.             continue_loop = going_to != coming_from || match_found;
  1586.             }
  1587.         } while (continue_loop && (col != curcol || row != currow));
  1588.  
  1589.         if (match_found <= 3)  /* DG */
  1590.         { /* no hole */
  1591.         color = bkcolor;
  1592.         reset_periodicity = 1;
  1593.         continue;
  1594.         }
  1595.  
  1596. /*
  1597. Fill in region by looping around again, filling lines to the left
  1598. whenever going_to is South or West
  1599. */
  1600.         trail_row = currow;
  1601.         trail_col = curcol;
  1602.         coming_from = West;
  1603.         going_to = East;
  1604.         do
  1605.         {
  1606.         match_found = FALSE;
  1607.         do
  1608.             {
  1609.             step_col_row();
  1610.             if (row >= currow
  1611.                 && col >= ixstart
  1612.                 && col <= ixstop
  1613.                 && row <= iystop
  1614.                 && getcolor(col,row) == trail_color)
  1615.                   /* getcolor() must be last */
  1616.             {
  1617.             if (going_to == South
  1618.                 || (going_to == West && coming_from != East))
  1619.                 { /* fill a row, but only once */
  1620.                 right = col;
  1621.                 while (--right >= ixstart && (color = getcolor(right,row)) == trail_color)
  1622.                 ; /* do nothing */
  1623.                 if (color == bkcolor) /* check last color */
  1624.                 {
  1625.                 left = right;
  1626.                 while (getcolor(--left,row) == bkcolor)
  1627.                       /* Should NOT be possible for left < ixstart */
  1628.                     ; /* do nothing */
  1629.                 left++; /* one pixel too far */
  1630.                 if (right == left) /* only one hole */
  1631.                     (*plot)(left,row,fillcolor_used);
  1632.                 else
  1633.                     { /* fill the line to the left */
  1634.                     length=right-left+1;
  1635.                     if (fillcolor_used != last_fillcolor_used || length > max_putline_length)
  1636.                     { /* only reset dstack if necessary */
  1637.                     memset(dstack,fillcolor_used,length);
  1638.                     last_fillcolor_used = fillcolor_used;
  1639.                     max_putline_length = length;
  1640.                     }
  1641.                     put_line(row,left,right,dstack);
  1642.                     /* here's where all the symmetry goes */
  1643.                     if (plot == putcolor)
  1644.                     kbdcount -= length >> 4; /* seems like a reasonable value */
  1645.                     else if (plot == symplot2) /* X-axis symmetry */
  1646.                     {
  1647.                     put_line(yystop-(row-yystart),left,right,dstack);
  1648.                     kbdcount -= length >> 3;
  1649.                     }
  1650.                     else if (plot == symplot2Y) /* Y-axis symmetry */
  1651.                     {
  1652.                     put_line(row,xxstop-(right-xxstart),xxstop-(left-xxstart),dstack);
  1653.                     kbdcount -= length >> 3;
  1654.                     }
  1655.                     else if (plot == symplot2J)  /* Origin symmetry */
  1656.                     {
  1657.                     put_line(yystop-(row-yystart),xxstop-(right-xxstart),xxstop-(left-xxstart),dstack);
  1658.                     kbdcount -= length >> 3;
  1659.                     }
  1660.                     else if (plot == symplot4) /* X-axis and Y-axis symmetry */
  1661.                     {
  1662.                     put_line(yystop-(row-yystart),left,right,dstack);
  1663.                     put_line(row,xxstop-(right-xxstart),xxstop-(left-xxstart),dstack);
  1664.                     put_line(yystop-(row-yystart),xxstop-(right-xxstart),xxstop-(left-xxstart),dstack);
  1665.                     kbdcount -= length >> 2;
  1666.                     }
  1667.                     else    /* cheap and easy way out */
  1668.                     {
  1669.                     int c;
  1670.                     for (c = left; c <= right; c++)  /* DG */
  1671.                         (*plot)(c,row,fillcolor_used);
  1672.                     kbdcount -= length >> 1;
  1673.                     }
  1674.                     }
  1675.                 } /* end of fill line */
  1676.  
  1677.                 if(--kbdcount<=0)
  1678.                 {
  1679.                 if(check_key())
  1680.                     {
  1681.                     if (iystop != yystop)
  1682.                        iystop = yystop - (currow - yystart); /* allow for sym */
  1683.                     add_worklist(xxstart,xxstop,currow,iystop,currow,0,worksym);
  1684.                     return(-1);
  1685.                     }
  1686.                 kbdcount=max_kbdcount;
  1687.                 }
  1688.                 }
  1689.             trail_row = row;
  1690.             trail_col = col;
  1691.             advance_match();
  1692.             match_found = TRUE;
  1693.             }
  1694.             else
  1695.             advance_no_match();
  1696.             } while (!match_found && going_to != coming_from);
  1697.  
  1698.         if (!match_found)
  1699.             { /* next one has to be a match */
  1700.             step_col_row();
  1701.             trail_row = row;
  1702.             trail_col = col;
  1703.             advance_match();
  1704.             }
  1705.         } while (trail_col != curcol || trail_row != currow);
  1706.         reset_periodicity = 1; /* reset after a trace/fill */
  1707.         color = bkcolor;
  1708.         }
  1709.     }
  1710.     return 0;
  1711.     }
  1712.  
  1713. /*******************************************************************/
  1714. /* take one step in the direction of going_to */
  1715. static void step_col_row()
  1716.     {
  1717.     switch (going_to)
  1718.     {
  1719.     case North:
  1720.         col = trail_col;
  1721.         row = trail_row - 1;
  1722.         break;
  1723.     case East:
  1724.         col = trail_col + 1;
  1725.         row = trail_row;
  1726.         break;
  1727.     case South:
  1728.         col = trail_col;
  1729.         row = trail_row + 1;
  1730.         break;
  1731.     case West:
  1732.         col = trail_col - 1;
  1733.         row = trail_row;
  1734.         break;
  1735.     }
  1736.     }
  1737.  
  1738. /******************* end of boundary trace method *******************/
  1739.  
  1740.  
  1741. /************************ super solid guessing *****************************/
  1742.  
  1743. /*
  1744.    I, Timothy Wegner, invented this solidguessing idea and implemented it in
  1745.    more or less the overall framework you see here.  I am adding this note
  1746.    now in a possibly vain attempt to secure my place in history, because
  1747.    Pieter Branderhorst has totally rewritten this routine, incorporating
  1748.    a *MUCH* more sophisticated algorithm.  His revised code is not only
  1749.    faster, but is also more accurate. Harrumph!
  1750. */
  1751.  
  1752. static int solidguess()
  1753. {
  1754.    int i,x,y,xlim,ylim,blocksize;
  1755.    unsigned int *pfxp0,*pfxp1;
  1756.    unsigned int u;
  1757.  
  1758.    guessplot=(plot!=putcolor && plot!=symplot2 && plot!=symplot2J);
  1759.    /* check if guessing at bottom & right edges is ok */
  1760.    bottom_guess = (plot == symplot2 || (plot == putcolor && iystop+1 == ydots));
  1761.    right_guess    = (plot == symplot2J
  1762.        || ((plot == putcolor || plot == symplot2) && ixstop+1 == xdots));
  1763.  
  1764.    i = maxblock = blocksize = ssg_blocksize();
  1765.    totpasses = 1;
  1766.    while ((i >>= 1) > 1) ++totpasses;
  1767.  
  1768.    /* ensure window top and left are on required boundary, treat window
  1769.      as larger than it really is if necessary (this is the reason symplot
  1770.      routines must check for > xdots/ydots before plotting sym points) */
  1771.    ixstart &= -1 - (maxblock-1);
  1772.    iystart = yybegin;
  1773.    iystart &= -1 - (maxblock-1);
  1774.  
  1775.    got_status = 1;
  1776.  
  1777.    if (workpass == 0) /* otherwise first pass already done */
  1778.    {
  1779.       /* first pass, calc every blocksize**2 pixel, quarter result & paint it */
  1780.       curpass = 1;
  1781.       if (iystart <= yystart) /* first time for this window, init it */
  1782.       {
  1783.      currow = 0;
  1784.      memset(&prefix[1][0][0],0,maxxblk*maxyblk*2); /* noskip flags off */
  1785.      reset_periodicity = 1;
  1786.      row=iystart;
  1787.      for(col=ixstart; col<=ixstop; col+=maxblock)
  1788.      { /* calc top row */
  1789.         if((*calctype)()== -1)
  1790.         {
  1791.            add_worklist(xxstart,xxstop,yystart,yystop,yybegin,0,worksym);
  1792.            goto exit_solidguess;
  1793.         }
  1794.         reset_periodicity = 0;
  1795.      }
  1796.       }
  1797.       else
  1798.      memset(&prefix[1][0][0],-1,maxxblk*maxyblk*2); /* noskip flags on */
  1799.       for(y=iystart; y<=iystop; y+=blocksize)
  1800.       {
  1801.      currow = y;
  1802.      i = 0;
  1803.      if(y+blocksize<=iystop)
  1804.      { /* calc the row below */
  1805.         row=y+blocksize;
  1806.         reset_periodicity = 1;
  1807.         for(col=ixstart; col<=ixstop; col+=maxblock)
  1808.         {
  1809.            if((i=(*calctype)()) == -1)
  1810.           break;
  1811.            reset_periodicity = 0;
  1812.         }
  1813.      }
  1814.      reset_periodicity = 1;
  1815.      if (i == -1 || guessrow(1,y,blocksize) != 0) /* interrupted? */
  1816.      {
  1817.         if (y < yystart)
  1818.            y = yystart;
  1819.         add_worklist(xxstart,xxstop,yystart,yystop,y,0,worksym);
  1820.         goto exit_solidguess;
  1821.      }
  1822.       }
  1823.  
  1824.       if (num_worklist) /* work list not empty, just do 1st pass */
  1825.       {
  1826.      add_worklist(xxstart,xxstop,yystart,yystop,yystart,1,worksym);
  1827.      goto exit_solidguess;
  1828.       }
  1829.       ++workpass;
  1830.       iystart = yystart & (-1 - (maxblock-1));
  1831.  
  1832.       /* calculate skip flags for skippable blocks */
  1833.       xlim=(ixstop+maxblock)/maxblock+1;
  1834.       ylim=((iystop+maxblock)/maxblock+15)/16+1;
  1835.       if(right_guess==0) /* no right edge guessing, zap border */
  1836.      for(y=0;y<=ylim;++y)
  1837.         prefix[1][y][xlim]= -1;
  1838.       if(bottom_guess==0) /* no bottom edge guessing, zap border */
  1839.       {
  1840.      i=(iystop+maxblock)/maxblock+1;
  1841.      y=i/16+1;
  1842.      i=1<<(i&15);
  1843.      for(x=0;x<=xlim;++x)
  1844.         prefix[1][y][x]|=i;
  1845.       }
  1846.       /* set each bit in prefix[0] to OR of it & surrounding 8 in prefix[1] */
  1847.       for(y=0;++y<ylim;)
  1848.       {
  1849.      pfxp0= &prefix[0][y][0];
  1850.      pfxp1= &prefix[1][y][0];
  1851.      for(x=0;++x<xlim;)
  1852.      {
  1853.         ++pfxp1;
  1854.         u= *(pfxp1-1)|*pfxp1|*(pfxp1+1);
  1855.         *(++pfxp0)=u|(u>>1)|(u<<1)
  1856.            |((*(pfxp1-(maxxblk+1))|*(pfxp1-maxxblk)|*(pfxp1-(maxxblk-1)))>>15)
  1857.           |((*(pfxp1+(maxxblk-1))|*(pfxp1+maxxblk)|*(pfxp1+(maxxblk+1)))<<15);
  1858.      }
  1859.       }
  1860.    }
  1861.    else /* first pass already done */
  1862.       memset(&prefix[0][0][0],-1,maxxblk*maxyblk*2); /* noskip flags on */
  1863.  
  1864.    /* remaining pass(es), halve blocksize & quarter each blocksize**2 */
  1865.    i = workpass;
  1866.    while (--i > 0) /* allow for already done passes */
  1867.       blocksize = blocksize>>1;
  1868.    reset_periodicity = 1;
  1869.    while((blocksize=blocksize>>1)>=2)
  1870.    {
  1871.       curpass = workpass + 1;
  1872.       for(y=iystart; y<=iystop; y+=blocksize)
  1873.       {
  1874.      currow = y;
  1875.      if(guessrow(0,y,blocksize)!=0)
  1876.      {
  1877.         if (y < yystart)
  1878.            y = yystart;
  1879.         add_worklist(xxstart,xxstop,yystart,yystop,y,workpass,worksym);
  1880.         goto exit_solidguess;
  1881.      }
  1882.       }
  1883.       ++workpass;
  1884.       if (num_worklist /* work list not empty, do one pass at a time */
  1885.       && blocksize>2) /* if 2, we just did last pass */
  1886.       {
  1887.      add_worklist(xxstart,xxstop,yystart,yystop,yystart,workpass,worksym);
  1888.      goto exit_solidguess;
  1889.       }
  1890.       iystart = yystart & (-1 - (maxblock-1));
  1891.    }
  1892.  
  1893.    exit_solidguess:
  1894.    return(0);
  1895. }
  1896.  
  1897. #define calcadot(c,x,y) { col=x; row=y; if((c=(*calctype)())== -1) return -1; }
  1898.  
  1899. static int _fastcall guessrow(int firstpass,int y,int blocksize)
  1900. {
  1901.    int x,i,j,color;
  1902.    int xplushalf,xplusblock;
  1903.    int ylessblock,ylesshalf,yplushalf,yplusblock;
  1904.    int       c21,c31,c41;     /* cxy is the color of pixel at (x,y) */
  1905.    int c12,c22,c32,c42;     /* where c22 is the topleft corner of */
  1906.    int c13,c23,c33;        /* the block being handled in current */
  1907.    int       c24,    c44;     /* iteration                  */
  1908.    int guessed23,guessed32,guessed33,guessed12,guessed13;
  1909.    int prev11,fix21,fix31;
  1910.    unsigned int *pfxptr,pfxmask;
  1911.  
  1912.    halfblock=blocksize>>1;
  1913.    i=y/maxblock;
  1914.    pfxptr= &prefix[firstpass][(i>>4)+1][ixstart/maxblock];
  1915.    pfxmask=1<<(i&15);
  1916.    ylesshalf=y-halfblock;
  1917.    ylessblock=y-blocksize; /* constants, for speed */
  1918.    yplushalf=y+halfblock;
  1919.    yplusblock=y+blocksize;
  1920.    prev11= -1;
  1921.    c24=c12=c13=c22=getcolor(ixstart,y);
  1922.    c31=c21=getcolor(ixstart,(y>0)?ylesshalf:0);
  1923.    if(yplusblock<=iystop)
  1924.       c24=getcolor(ixstart,yplusblock);
  1925.    else if(bottom_guess==0)
  1926.       c24= -1;
  1927.    guessed12=guessed13=0;
  1928.  
  1929.    for(x=ixstart; x<=ixstop;)  /* increment at end, or when doing continue */
  1930.    {
  1931.       if((x&(maxblock-1))==0)  /* time for skip flag stuff */
  1932.       {
  1933.      ++pfxptr;
  1934.      if(firstpass==0 && (*pfxptr&pfxmask)==0)  /* check for fast skip */
  1935.      {
  1936.         /* next useful in testing to make skips visible */
  1937.         /*
  1938.                if(halfblock==1)
  1939.                {
  1940.                   (*plot)(x+1,y,0); (*plot)(x,y+1,0); (*plot)(x+1,y+1,0);
  1941.                   }
  1942.              */
  1943.         x+=maxblock;
  1944.         prev11=c31=c21=c24=c12=c13=c22;
  1945.         guessed12=guessed13=0;
  1946.         continue;
  1947.      }
  1948.       }
  1949.  
  1950.       if(firstpass)  /* 1st pass, paint topleft corner */
  1951.      plotblock(0,x,y,c22);
  1952.       /* setup variables */
  1953.       xplusblock=(xplushalf=x+halfblock)+halfblock;
  1954.       if(xplushalf>ixstop)
  1955.       {
  1956.      if(right_guess==0)
  1957.         c31= -1;
  1958.       }
  1959.       else if(y>0)
  1960.      c31=getcolor(xplushalf,ylesshalf);
  1961.       if(xplusblock<=ixstop)
  1962.       {
  1963.      if(yplusblock<=iystop)
  1964.         c44=getcolor(xplusblock,yplusblock);
  1965.      c41=getcolor(xplusblock,(y>0)?ylesshalf:0);
  1966.      c42=getcolor(xplusblock,y);
  1967.       }
  1968.       else if(right_guess==0)
  1969.      c41=c42=c44= -1;
  1970.       if(yplusblock>iystop)
  1971.      c44=(bottom_guess)?c42:-1;
  1972.  
  1973.       /* guess or calc the remaining 3 quarters of current block */
  1974.       guessed23=guessed32=guessed33=1;
  1975.       c23=c32=c33=c22;
  1976.       if(yplushalf>iystop)
  1977.       {
  1978.      if(bottom_guess==0)
  1979.         c23=c33= -1;
  1980.      guessed23=guessed33= -1;
  1981.       }
  1982.       if(xplushalf>ixstop)
  1983.       {
  1984.      if(right_guess==0)
  1985.         c32=c33= -1;
  1986.      guessed32=guessed33= -1;
  1987.       }
  1988.       while(1) /* go around till none of 23,32,33 change anymore */
  1989.       {
  1990.      if(guessed33>0
  1991.          && (c33!=c44 || c33!=c42 || c33!=c24 || c33!=c32 || c33!=c23))
  1992.      {
  1993.         calcadot(c33,xplushalf,yplushalf);
  1994.         guessed33=0;
  1995.      }
  1996.      if(guessed32>0
  1997.          && (c32!=c33 || c32!=c42 || c32!=c31 || c32!=c21
  1998.          || c32!=c41 || c32!=c23))
  1999.      {
  2000.         calcadot(c32,xplushalf,y);
  2001.         guessed32=0;
  2002.         continue;
  2003.      }
  2004.      if(guessed23>0
  2005.          && (c23!=c33 || c23!=c24 || c23!=c13 || c23!=c12 || c23!=c32))
  2006.      {
  2007.         calcadot(c23,x,yplushalf);
  2008.         guessed23=0;
  2009.         continue;
  2010.      }
  2011.      break;
  2012.       }
  2013.  
  2014.       if(firstpass) /* note whether any of block's contents were calculated */
  2015.      if(guessed23==0 || guessed32==0 || guessed33==0)
  2016.         *pfxptr|=pfxmask;
  2017.  
  2018.       if(halfblock>1) /* not last pass, check if something to display */
  2019.      if(firstpass)    /* display guessed corners, fill in block */
  2020.      {
  2021.         if(guessplot)
  2022.         {
  2023.            if(guessed23>0)
  2024.           (*plot)(x,yplushalf,c23);
  2025.            if(guessed32>0)
  2026.           (*plot)(xplushalf,y,c32);
  2027.            if(guessed33>0)
  2028.           (*plot)(xplushalf,yplushalf,c33);
  2029.         }
  2030.         plotblock(1,x,yplushalf,c23);
  2031.         plotblock(0,xplushalf,y,c32);
  2032.         plotblock(1,xplushalf,yplushalf,c33);
  2033.      }
  2034.      else  /* repaint changed blocks */
  2035.      {
  2036.         if(c23!=c22)
  2037.            plotblock(-1,x,yplushalf,c23);
  2038.         if(c32!=c22)
  2039.            plotblock(-1,xplushalf,y,c32);
  2040.         if(c33!=c22)
  2041.            plotblock(-1,xplushalf,yplushalf,c33);
  2042.      }
  2043.  
  2044.       /* check if some calcs in this block mean earlier guesses need fixing */
  2045.       fix21=((c22!=c12 || c22!=c32)
  2046.       && c21==c22 && c21==c31 && c21==prev11
  2047.       && y>0
  2048.       && (x==ixstart || c21==getcolor(x-halfblock,ylessblock))
  2049.       && (xplushalf>ixstop || c21==getcolor(xplushalf,ylessblock))
  2050.       && c21==getcolor(x,ylessblock));
  2051.       fix31=(c22!=c32
  2052.       && c31==c22 && c31==c42 && c31==c21 && c31==c41
  2053.       && y>0 && xplushalf<=ixstop
  2054.       && c31==getcolor(xplushalf,ylessblock)
  2055.       && (xplusblock>ixstop || c31==getcolor(xplusblock,ylessblock))
  2056.       && c31==getcolor(x,ylessblock));
  2057.       prev11=c31; /* for next time around */
  2058.       if(fix21)
  2059.       {
  2060.      calcadot(c21,x,ylesshalf);
  2061.      if(halfblock>1 && c21!=c22)
  2062.         plotblock(-1,x,ylesshalf,c21);
  2063.       }
  2064.       if(fix31)
  2065.       {
  2066.      calcadot(c31,xplushalf,ylesshalf);
  2067.      if(halfblock>1 && c31!=c22)
  2068.         plotblock(-1,xplushalf,ylesshalf,c31);
  2069.       }
  2070.       if(c23!=c22)
  2071.       {
  2072.      if(guessed12)
  2073.      {
  2074.         calcadot(c12,x-halfblock,y);
  2075.         if(halfblock>1 && c12!=c22)
  2076.            plotblock(-1,x-halfblock,y,c12);
  2077.      }
  2078.      if(guessed13)
  2079.      {
  2080.         calcadot(c13,x-halfblock,yplushalf);
  2081.         if(halfblock>1 && c13!=c22)
  2082.            plotblock(-1,x-halfblock,yplushalf,c13);
  2083.      }
  2084.       }
  2085.       c22=c42;
  2086.       c24=c44;
  2087.       c13=c33;
  2088.       c31=c21=c41;
  2089.       c12=c32;
  2090.       guessed12=guessed32;
  2091.       guessed13=guessed33;
  2092.       x+=blocksize;
  2093.    } /* end x loop */
  2094.  
  2095.    if(firstpass==0 || guessplot) return 0;
  2096.  
  2097.    /* paint rows the fast way */
  2098.    for(i=0;i<halfblock;++i)
  2099.    {
  2100.       if((j=y+i)<=iystop)
  2101.      put_line(j,xxstart,ixstop,&dstack[xxstart]);
  2102.       if((j=y+i+halfblock)<=iystop)
  2103.      put_line(j,xxstart,ixstop,&dstack[xxstart+MAXPIXELS]);
  2104.       if(keypressed()) return -1;
  2105.    }
  2106.    if(plot!=putcolor)  /* symmetry, just vertical & origin the fast way */
  2107.    {
  2108.       if(plot==symplot2J) /* origin sym, reverse lines */
  2109.      for(i=(ixstop+xxstart+1)/2;--i>=xxstart;)
  2110.      {
  2111.         color=dstack[i];
  2112.         dstack[i]=dstack[j=ixstop-(i-xxstart)];
  2113.         dstack[j]=color;
  2114.         j+=MAXPIXELS;
  2115.         color=dstack[i+MAXPIXELS];
  2116.         dstack[i+MAXPIXELS]=dstack[j];
  2117.         dstack[j]=color;
  2118.      }
  2119.       for(i=0;i<halfblock;++i)
  2120.       {
  2121.      if((j=yystop-(y+i-yystart))>iystop && j<ydots)
  2122.         put_line(j,xxstart,ixstop,&dstack[xxstart]);
  2123.      if((j=yystop-(y+i+halfblock-yystart))>iystop && j<ydots)
  2124.         put_line(j,xxstart,ixstop,&dstack[xxstart+MAXPIXELS]);
  2125.      if(keypressed()) return -1;
  2126.       }
  2127.    }
  2128.    return 0;
  2129. }
  2130.  
  2131. static void _fastcall plotblock(int buildrow,int x,int y,int color)
  2132. {
  2133.    int i,xlim,ylim;
  2134.    if((xlim=x+halfblock)>ixstop)
  2135.       xlim=ixstop+1;
  2136.    if(buildrow>=0 && guessplot==0) /* save it for later put_line */
  2137.    {
  2138.       if(buildrow==0)
  2139.      for(i=x;i<xlim;++i)
  2140.         dstack[i]=color;
  2141.       else
  2142.      for(i=x;i<xlim;++i)
  2143.         dstack[i+MAXPIXELS]=color;
  2144.       if (x>=xxstart) /* when x reduced for alignment, paint those dots too */
  2145.      return; /* the usual case */
  2146.    }
  2147.    /* paint it */
  2148.    if((ylim=y+halfblock)>iystop)
  2149.    {
  2150.       if(y>iystop)
  2151.      return;
  2152.       ylim=iystop+1;
  2153.    }
  2154.    for(i=x;++i<xlim;)
  2155.       (*plot)(i,y,color); /* skip 1st dot on 1st row */
  2156.    while(++y<ylim)
  2157.       for(i=x;i<xlim;++i)
  2158.      (*plot)(i,y,color);
  2159. }
  2160.  
  2161.  
  2162. /************************* symmetry plot setup ************************/
  2163.  
  2164. static int _fastcall xsym_split(int xaxis_row,int xaxis_between)
  2165. {
  2166.    int i;
  2167.    if ((worksym&0x11) == 0x10) /* already decided not sym */
  2168.       return(1);
  2169.    if ((worksym&1) != 0) /* already decided on sym */
  2170.       iystop = (yystart+yystop)/2;
  2171.    else /* new window, decide */
  2172.    {
  2173.       worksym |= 0x10;
  2174.       if (xaxis_row <= yystart || xaxis_row >= yystop)
  2175.      return(1); /* axis not in window */
  2176.       i = xaxis_row + (xaxis_row - yystart);
  2177.       if (xaxis_between)
  2178.      ++i;
  2179.       if (i > yystop) /* split into 2 pieces, bottom has the symmetry */
  2180.       {
  2181.      if (num_worklist >= MAXCALCWORK-1) /* no room to split */
  2182.         return(1);
  2183.      iystop = xaxis_row - (yystop - xaxis_row);
  2184.      if (!xaxis_between)
  2185.         --iystop;
  2186.      add_worklist(xxstart,xxstop,iystop+1,yystop,iystop+1,workpass,0);
  2187.      yystop = iystop;
  2188.      return(1); /* tell set_symmetry no sym for current window */
  2189.       }
  2190.       if (i < yystop) /* split into 2 pieces, top has the symmetry */
  2191.       {
  2192.      if (num_worklist >= MAXCALCWORK-1) /* no room to split */
  2193.         return(1);
  2194.      add_worklist(xxstart,xxstop,i+1,yystop,i+1,workpass,0);
  2195.      yystop = i;
  2196.       }
  2197.       iystop = xaxis_row;
  2198.       worksym |= 1;
  2199.    }
  2200.    symmetry = 0;
  2201.    return(0); /* tell set_symmetry its a go */
  2202. }
  2203.  
  2204. static int _fastcall ysym_split(int yaxis_col,int yaxis_between)
  2205. {
  2206.    int i;
  2207.    if ((worksym&0x22) == 0x20) /* already decided not sym */
  2208.       return(1);
  2209.    if ((worksym&2) != 0) /* already decided on sym */
  2210.       ixstop = (xxstart+xxstop)/2;
  2211.    else /* new window, decide */
  2212.    {
  2213.       worksym |= 0x20;
  2214.       if (yaxis_col <= xxstart || yaxis_col >= xxstop)
  2215.      return(1); /* axis not in window */
  2216.       i = yaxis_col + (yaxis_col - xxstart);
  2217.       if (yaxis_between)
  2218.      ++i;
  2219.       if (i > xxstop) /* split into 2 pieces, right has the symmetry */
  2220.       {
  2221.      if (num_worklist >= MAXCALCWORK-1) /* no room to split */
  2222.         return(1);
  2223.      ixstop = yaxis_col - (xxstop - yaxis_col);
  2224.      if (!yaxis_between)
  2225.         --ixstop;
  2226.      add_worklist(ixstop+1,xxstop,yystart,yystop,yystart,workpass,0);
  2227.      xxstop = ixstop;
  2228.      return(1); /* tell set_symmetry no sym for current window */
  2229.       }
  2230.       if (i < xxstop) /* split into 2 pieces, left has the symmetry */
  2231.       {
  2232.      if (num_worklist >= MAXCALCWORK-1) /* no room to split */
  2233.         return(1);
  2234.      add_worklist(i+1,xxstop,yystart,yystop,yystart,workpass,0);
  2235.      xxstop = i;
  2236.       }
  2237.       ixstop = yaxis_col;
  2238.       worksym |= 2;
  2239.    }
  2240.    symmetry = 0;
  2241.    return(0); /* tell set_symmetry its a go */
  2242. }
  2243.  
  2244. static void _fastcall setsymmetry(int sym, int uselist) /* set up proper symmetrical plot functions */
  2245. {
  2246.    extern int forcesymmetry;
  2247.    int i;
  2248.    int parmszero;
  2249.    int xaxis_row, yaxis_col;         /* pixel number for origin */
  2250.    int xaxis_between, yaxis_between; /* if axis between 2 pixels, not on one */
  2251.    double ftemp;
  2252.    symmetry = 1;
  2253.    TranspSymmetry = sym;         /* for tp3d.c, MCP 6-1-90 */
  2254.    if(sym == NOPLOT && forcesymmetry == 999)
  2255.    {
  2256.       plot = noplot;
  2257.       return;
  2258.    }
  2259.    /* NOTE: 16-bit potential disables symmetry */
  2260.    /* also any decomp= option and any inversion not about the origin */
  2261.    /* also any rotation other than 180deg and any off-axis stretch */
  2262.    if ((potflag && pot16bit) || (invert && inversion[2] != 0.0)
  2263.        || decomp[0] != 0
  2264.        || xxmin!=xx3rd || yymin!=yy3rd)
  2265.       return;
  2266.    if(sym != XAXIS && sym != XAXIS_NOPARM && inversion[1] != 0.0 && forcesymmetry == 999)
  2267.       return;
  2268.    if(forcesymmetry != 999)
  2269.       sym = forcesymmetry;
  2270.    parmszero = (parm.x == 0.0 && parm.y == 0.0 && useinitorbit != 1);
  2271.    xaxis_row = yaxis_col = -1;
  2272.    if (fabs(yymin+yymax) < fabs(yymin)+fabs(yymax)) /* axis is on screen */
  2273.    {
  2274.       ftemp = (0.0-yymax) / (yymin-yymax) * (ydots-1) + 0.25;
  2275.       xaxis_row = ftemp;
  2276.       xaxis_between = (ftemp - xaxis_row >= 0.5);
  2277.       if (uselist == 0 && (!xaxis_between || (xaxis_row+1)*2 != ydots))
  2278.      xaxis_row = -1; /* can't split screen, so dead center or not at all */
  2279.    }
  2280.    if (fabs(xxmin+xxmax) < fabs(xxmin)+fabs(xxmax)) /* axis is on screen */
  2281.    {
  2282.       ftemp = (0.0-xxmin) / (xxmax-xxmin) * (xdots-1) + 0.25;
  2283.       yaxis_col = ftemp;
  2284.       yaxis_between = (ftemp - yaxis_col >= 0.5);
  2285.       if (uselist == 0 && (!yaxis_between || (yaxis_col+1)*2 != xdots))
  2286.      yaxis_col = -1; /* can't split screen, so dead center or not at all */
  2287.    }
  2288.    switch(sym)         /* symmetry switch */
  2289.    {
  2290.    case XAXIS_NOREAL:     /* X-axis Symmetry (no real param) */
  2291.       if (parm.x != 0.0) break;
  2292.       goto xsym;
  2293.    case XAXIS_NOIMAG:     /* X-axis Symmetry (no imag param) */
  2294.       if (parm.y != 0.0) break;
  2295.       goto xsym;
  2296.    case XAXIS_NOPARM:                 /* X-axis Symmetry  (no params)*/
  2297.       if (!parmszero)
  2298.      break;
  2299.       xsym:
  2300.    case XAXIS:                 /* X-axis Symmetry */
  2301.       if (xsym_split(xaxis_row,xaxis_between) == 0)
  2302.      if(basin)
  2303.         plot = symplot2basin;
  2304.      else
  2305.         plot = symplot2;
  2306.       break;
  2307.    case YAXIS_NOPARM:                 /* Y-axis Symmetry (No Parms)*/
  2308.       if (!parmszero)
  2309.      break;
  2310.    case YAXIS:                 /* Y-axis Symmetry */
  2311.       if (ysym_split(yaxis_col,yaxis_between) == 0)
  2312.      plot = symplot2Y;
  2313.       break;
  2314.    case XYAXIS_NOPARM:                 /* X-axis AND Y-axis Symmetry (no parms)*/
  2315.       if(!parmszero)
  2316.      break;
  2317.    case XYAXIS:              /* X-axis AND Y-axis Symmetry */
  2318.       xsym_split(xaxis_row,xaxis_between);
  2319.       ysym_split(yaxis_col,yaxis_between);
  2320.       switch (worksym & 3)
  2321.       {
  2322.       case 1: /* just xaxis symmetry */
  2323.      if(basin)
  2324.         plot = symplot2basin;
  2325.      else
  2326.         plot = symplot2 ;
  2327.      break;
  2328.       case 2: /* just yaxis symmetry */
  2329.      if (basin) /* got no routine for this case */
  2330.      {
  2331.         ixstop = xxstop; /* fix what split should not have done */
  2332.         symmetry = 1;
  2333.      }
  2334.      else
  2335.         plot = symplot2Y;
  2336.      break;
  2337.       case 3: /* both axes */
  2338.      if(basin)
  2339.         plot = symplot4basin;
  2340.      else
  2341.         plot = symplot4 ;
  2342.       }
  2343.       break;
  2344.    case ORIGIN_NOPARM:                 /* Origin Symmetry (no parms)*/
  2345.       if (!parmszero)
  2346.      break;
  2347.    case ORIGIN:              /* Origin Symmetry */
  2348.       originsym:
  2349.       if ( xsym_split(xaxis_row,xaxis_between) == 0
  2350.       && ysym_split(yaxis_col,yaxis_between) == 0)
  2351.       {
  2352.      plot = symplot2J;
  2353.      ixstop = xxstop; /* didn't want this changed */
  2354.       }
  2355.       else
  2356.       {
  2357.      iystop = yystop; /* in case first split worked */
  2358.      symmetry = 1;
  2359.      worksym = 0x30; /* let it recombine with others like it */
  2360.       }
  2361.       break;
  2362.    case PI_SYM_NOPARM:
  2363.       if (!parmszero)
  2364.      break;
  2365.    case PI_SYM:              /* PI symmetry */
  2366.       if(invert && forcesymmetry == 999)
  2367.      goto originsym;
  2368.       plot = symPIplot ;
  2369.       symmetry = 0;
  2370.       if ( xsym_split(xaxis_row,xaxis_between) == 0
  2371.       && ysym_split(yaxis_col,yaxis_between) == 0)
  2372.      if(parm.y == 0.0)
  2373.         plot = symPIplot4J; /* both axes */
  2374.      else
  2375.         plot = symPIplot2J; /* origin */
  2376.       else
  2377.       {
  2378.      iystop = yystop; /* in case first split worked */
  2379.      worksym = 0x30;  /* don't mark pisym as ysym, just do it unmarked */
  2380.       }
  2381.       pixelpi = (PI/fabs(xxmax-xxmin))*xdots; /* PI in pixels */
  2382.       if ( (ixstop = xxstart+pixelpi-1 ) > xxstop)
  2383.      ixstop = xxstop;
  2384.       if (plot == symPIplot4J && ixstop > (i = (xxstart+xxstop)/2))
  2385.      ixstop = i;
  2386.       break;
  2387.    default:             /* no symmetry */
  2388.       break;
  2389.    }
  2390. }
  2391.  
  2392.  
  2393. /***************** standalone engine for "test" ********************/
  2394.  
  2395. int test()
  2396. {
  2397.    int startrow,startpass,numpasses;
  2398.    startrow = startpass = 0;
  2399.    if (resuming)
  2400.    {
  2401.       start_resume();
  2402.       get_resume(sizeof(int),&startrow,sizeof(int),&startpass,0);
  2403.       end_resume();
  2404.    }
  2405.    if(teststart()) /* assume it was stand-alone, doesn't want passes logic */
  2406.       return(0);
  2407.    numpasses = (stdcalcmode == '1') ? 0 : 1;
  2408.    for (passes=startpass; passes <= numpasses ; passes++)
  2409.    {
  2410.       for (row = startrow; row <= iystop; row=row+1+numpasses)
  2411.       {
  2412.      register int col;
  2413.      for (col = 0; col <= ixstop; col++)       /* look at each point on screen */
  2414.      {
  2415.         register color;
  2416.         init.y = dy0[row]+dy1[col];
  2417.         init.x = dx0[col]+dx1[row];
  2418.         if(check_key())
  2419.         {
  2420.            testend();
  2421.            alloc_resume(20,1);
  2422.            put_resume(sizeof(int),&row,sizeof(int),&passes,0);
  2423.            return(-1);
  2424.         }
  2425.         color = testpt(init.x,init.y,parm.x,parm.y,maxit,inside);
  2426.         if (color >= colors) /* avoid trouble if color is 0 */
  2427.            if (colors < 16)
  2428.           color &= andcolor;
  2429.            else
  2430.           color = ((color-1) % andcolor) + 1; /* skip color zero */
  2431.         (*plot)(col,row,color);
  2432.         if(numpasses && (passes == 0))
  2433.            (*plot)(col,row+1,color);
  2434.      }
  2435.       }
  2436.       startrow = passes + 1;
  2437.    }
  2438.    testend();
  2439.    return(0);
  2440. }
  2441.  
  2442.  
  2443. /***************** standalone engine for "plasma" ********************/
  2444.  
  2445. static int iparmx;                /* iparmx = parm.x * 16 */
  2446. static int shiftvalue;                /* shift based on #colors */
  2447. extern int max_colors;
  2448. static int recur1=1;
  2449. static int pcolors;
  2450.  
  2451. typedef struct palett
  2452. {
  2453.    BYTE red;
  2454.    BYTE green;
  2455.    BYTE blue;
  2456. }
  2457. Palettetype;
  2458.  
  2459. extern Palettetype dacbox[256];
  2460. static int plasma_check;            /* to limit kbd checking */
  2461.  
  2462. static int _fastcall adjust(int xa,int ya,int x,int y,int xb,int yb)
  2463. {
  2464.    long pseudorandom;
  2465.    pseudorandom = ((long)iparmx)*((rand15()-16383));
  2466.    pseudorandom = pseudorandom*(abs(xa-xb)+abs(ya-yb));
  2467.    pseudorandom = pseudorandom >> shiftvalue;
  2468.    pseudorandom = ((getcolor(xa,ya)+getcolor(xb,yb)+1)>>1)+pseudorandom;
  2469.    if (pseudorandom <    1) pseudorandom =   1;
  2470.    if (pseudorandom >= pcolors) pseudorandom = pcolors-1;
  2471.    putcolor(x,y,(int)pseudorandom);
  2472.    return((int)pseudorandom);
  2473. }
  2474.  
  2475.  
  2476. static int _fastcall new_adjust(int xa,int ya,int x,int y,int xb,int yb)
  2477. {
  2478.    long pseudorandom;
  2479.    pseudorandom = ((long)iparmx)*((rand15()-16383));
  2480.    pseudorandom = pseudorandom * recur1;
  2481.    pseudorandom = pseudorandom >> shiftvalue;
  2482.    pseudorandom = ((getcolor(xa,ya)+getcolor(xb,yb)+1)>>1)+pseudorandom;
  2483.    if (pseudorandom <    1) pseudorandom =   1;
  2484.    if (pseudorandom >= pcolors) pseudorandom = pcolors-1;
  2485.    putcolor(x,y,(int)pseudorandom);
  2486.    return((int)pseudorandom);
  2487. }
  2488.  
  2489.  
  2490. static int _fastcall new_subD (int x1,int y1,int x2,int y2, int recur)
  2491. {
  2492.  
  2493.    int x,y;
  2494.    int v;
  2495.    int nx1;
  2496.    int nx;
  2497.    int ny1, ny;
  2498.    int i;
  2499.  
  2500.     struct sub {
  2501.     BYTE t; /* top of stack */
  2502.     int v[16]; /* subdivided value */
  2503.     BYTE r[16];  /* recursion level */
  2504.     };
  2505.  
  2506.     static struct sub subx, suby;
  2507.  
  2508.     recur1=1;
  2509.     for (i=1;i<=recur;i++)
  2510.     recur1 = recur1 * 2;
  2511.     recur1=320/recur1;
  2512.     suby.t = 2;
  2513.     ny    = suby.v[0] = y2;
  2514.     ny1 = suby.v[2] = y1;
  2515.     suby.r[0] = suby.r[2] = 0;
  2516.     suby.r[1] = 1;
  2517.     y = suby.v[1] = (ny1 + ny) >> 1;
  2518.  
  2519.     while (suby.t >= 1)
  2520.     {
  2521.     if ((++plasma_check & 0x0f) == 1)
  2522.         if(check_key())
  2523.         {
  2524. /*   naah, we don't want to flush this key!!!
  2525.         getch();
  2526. */
  2527.         plasma_check--;
  2528.         return(1);
  2529.         }
  2530.     while (suby.r[suby.t-1] < recur)
  2531.     {
  2532.         /*    1.  Create new entry at top of the stack  */
  2533.  
  2534.         /*    2.  Copy old top value to new top value.
  2535.             This is largest y value.          */
  2536.  
  2537.         /*    3.  Smallest y is now old mid point      */
  2538.  
  2539.         /*    4.  Set new mid point recursion level      */
  2540.  
  2541.         /*    5.  New mid point value is average
  2542.             of largest and smallest          */
  2543.  
  2544.         suby.t++;
  2545.         ny1  = suby.v[suby.t] = suby.v[suby.t-1];
  2546.         ny     = suby.v[suby.t-2];
  2547.         suby.r[suby.t] = suby.r[suby.t-1];
  2548.         y     = suby.v[suby.t-1]   = (ny1 + ny) >> 1;
  2549.         suby.r[suby.t-1]   = (int)max(suby.r[suby.t],
  2550.                         suby.r[suby.t-2])+1;
  2551.     }
  2552.     subx.t = 2;
  2553.     nx  = subx.v[0] = x2;
  2554.     nx1 = subx.v[2] = x1;
  2555.     subx.r[0] = subx.r[2] = 0;
  2556.     subx.r[1] = 1;
  2557.     x = subx.v[1] = (nx1 + nx) >> 1;
  2558.  
  2559.     while (subx.t >= 1)
  2560.     {
  2561.         while (subx.r[subx.t-1] < recur)
  2562.         {
  2563.         subx.t++; /* move the top ofthe stack up 1 */
  2564.         nx1  = subx.v[subx.t] = subx.v[subx.t-1];
  2565.         nx   = subx.v[subx.t-2];
  2566.         subx.r[subx.t] = subx.r[subx.t-1];
  2567.         x    = subx.v[subx.t-1]   = (nx1 + nx) >> 1;
  2568.         subx.r[subx.t-1]   = (int)max(subx.r[subx.t],
  2569.                         subx.r[subx.t-2])+1;
  2570.  
  2571.     }
  2572.  
  2573.        if ((i = getcolor(nx, y)) == 0)
  2574.            i = new_adjust(nx,ny1,nx,y ,nx,ny);
  2575.            v = i;
  2576.        if ((i = getcolor(x, ny)) == 0)
  2577.            i = new_adjust(nx1,ny,x ,ny,nx,ny);
  2578.            v += i;
  2579.        if(getcolor(x,y) == 0)
  2580.        {
  2581.            if ((i = getcolor(x, ny1)) == 0)
  2582.            i = new_adjust(nx1,ny1,x ,ny1,nx,ny1);
  2583.            v += i;
  2584.            if ((i = getcolor(nx1, y)) == 0)
  2585.            i = new_adjust(nx1,ny1,nx1,y ,nx1,ny);
  2586.            v += i;
  2587.            putcolor(x,y,(v + 2) >> 2);
  2588.        }
  2589.  
  2590.        if (subx.r[subx.t-1] == recur) subx.t = subx.t - 2;
  2591.     }
  2592.  
  2593.     if (suby.r[suby.t-1] == recur) suby.t = suby.t - 2;
  2594.     }
  2595.     return(0);
  2596. }
  2597.  
  2598.  
  2599. static void _fastcall subDivide(int x1,int y1,int x2,int y2)
  2600. {
  2601.    int x,y;
  2602.    int v,i;
  2603.  
  2604.    if ((++plasma_check & 0x7f) == 1)
  2605.       if(check_key())
  2606.       {
  2607.      plasma_check--;
  2608.      return;
  2609.       }
  2610.    if(x2-x1<2 && y2-y1<2)
  2611.       return;
  2612.  
  2613.    x = (x1+x2)>>1;
  2614.    y = (y1+y2)>>1;
  2615.    if((v=getcolor(x,y1)) == 0)
  2616.        v=adjust(x1,y1,x ,y1,x2,y1);
  2617.    i=v;
  2618.    if((v=getcolor(x2,y)) == 0)
  2619.        v=adjust(x2,y1,x2,y ,x2,y2);
  2620.    i+=v;
  2621.    if((v=getcolor(x,y2)) == 0)
  2622.        v=adjust(x1,y2,x ,y2,x2,y2);
  2623.    i+=v;
  2624.    if((v=getcolor(x1,y)) == 0)
  2625.        v=adjust(x1,y1,x1,y ,x1,y2);
  2626.    i+=v;
  2627.  
  2628.    if(getcolor(x,y) == 0)
  2629.       putcolor(x,y,(i+2)>>2);
  2630.  
  2631.    subDivide(x1,y1,x ,y);
  2632.    subDivide(x ,y1,x2,y);
  2633.    subDivide(x ,y ,x2,y2);
  2634.    subDivide(x1,y ,x ,y2);
  2635. }
  2636.  
  2637. int plasma()
  2638. {
  2639. int i,k;
  2640.  
  2641.    plasma_check = 0;
  2642.  
  2643.    if(colors < 4) {
  2644. static char far plasmamsg[]={"\
  2645. Plasma Clouds can currently only be run in a 4-or-more-color video\n\
  2646. mode (and color-cycled only on VGA adapters [or EGA adapters in their\n\
  2647. 640x350x16 mode])."};
  2648.       stopmsg(0,plasmamsg);
  2649.       return(-1);
  2650.    }
  2651.    iparmx = param[0] * 8;
  2652.    if (parm.x <= 0.0) iparmx = 16;
  2653.    if (parm.x >= 100) iparmx = 800;
  2654.  
  2655.    if ((!rflag) && param[2] == 1)
  2656.        --rseed;
  2657.    if (param[2] != 0 && param[2] != 1)
  2658.       rseed = param[2];
  2659.  
  2660.    srand(rseed);
  2661.    if (!rflag) ++rseed;
  2662.  
  2663.    if (colors == 256)            /* set the (256-color) palette */
  2664.       set_Plasma_palette();        /* skip this if < 256 colors */
  2665.  
  2666.    if (colors > 16)            /* adjust based on # of colors */
  2667.       shiftvalue = 18;
  2668.    else {
  2669.       if (colors > 4)
  2670.      shiftvalue = 22;
  2671.       else {
  2672.      if (colors > 2)
  2673.         shiftvalue = 24;
  2674.      else
  2675.         shiftvalue = 25;
  2676.       }
  2677.    }
  2678.  
  2679.    pcolors = min(colors, max_colors);
  2680.    putcolor(      0,      0,1+(((rand15()/pcolors)*(pcolors-1))>>(shiftvalue-11)));
  2681.    putcolor(xdots-1,      0,1+(((rand15()/pcolors)*(pcolors-1))>>(shiftvalue-11)));
  2682.    putcolor(xdots-1,ydots-1,1+(((rand15()/pcolors)*(pcolors-1))>>(shiftvalue-11)));
  2683.    putcolor(      0,ydots-1,1+(((rand15()/pcolors)*(pcolors-1))>>(shiftvalue-11)));
  2684.  
  2685.    if (param[1] == 0)
  2686.    {
  2687.       subDivide(0,0,xdots-1,ydots-1);
  2688.    }
  2689.    else
  2690.    {
  2691.       recur1 = i = k = 1;
  2692.       while(new_subD(0,0,xdots-1,ydots-1,i)==0)
  2693.       {
  2694.      k = k * 2;
  2695.      if (k    >(int)max(xdots-1,ydots-1))
  2696.          break;
  2697.      if (check_key())
  2698.          return(1);
  2699.      i++;
  2700.       }
  2701.    }
  2702.  
  2703.    if (! check_key())
  2704.       return(0);
  2705.    else
  2706.       return(1);
  2707. }
  2708.  
  2709. static void set_Plasma_palette()
  2710. {
  2711.    extern char far *mapdacbox;
  2712.    static Palettetype Red    = {
  2713.       63, 0, 0     };
  2714.    static Palettetype Green  = {
  2715.       0,63, 0    };
  2716.    static Palettetype Blue   = {
  2717.       0, 0,63    };
  2718.    int i;
  2719.  
  2720.    if (mapdacbox) return;        /* map= specified */
  2721.  
  2722.    dacbox[0].red  = 0 ;
  2723.    dacbox[0].green= 0 ;
  2724.    dacbox[0].blue = 0 ;
  2725.    for(i=1;i<=85;i++)
  2726.    {
  2727.       dacbox[i].red   = (i*Green.red   + (86-i)*Blue.red)/85;
  2728.       dacbox[i].green = (i*Green.green + (86-i)*Blue.green)/85;
  2729.       dacbox[i].blue  = (i*Green.blue  + (86-i)*Blue.blue)/85;
  2730.  
  2731.       dacbox[i+85].red     = (i*Red.red    + (86-i)*Green.red)/85;
  2732.       dacbox[i+85].green = (i*Red.green + (86-i)*Green.green)/85;
  2733.       dacbox[i+85].blue  = (i*Red.blue    + (86-i)*Green.blue)/85;
  2734.  
  2735.       dacbox[i+170].red   = (i*Blue.red   + (86-i)*Red.red)/85;
  2736.       dacbox[i+170].green = (i*Blue.green + (86-i)*Red.green)/85;
  2737.       dacbox[i+170].blue  = (i*Blue.blue  + (86-i)*Red.blue)/85;
  2738.    }
  2739.    SetTgaColors();    /* TARGA 3 June 89  j mclain */
  2740.    spindac(0,1);
  2741. }
  2742.  
  2743. /***************** standalone engine for "diffusion" ********************/
  2744.  
  2745. #define RANDOM(x)  (rand()%(x))
  2746.  
  2747. /* This constant assumes that rand() returns a value from 0 to 32676 */
  2748. #define FOURPI  (long)(4*PI*(double)(1L << 16))
  2749.  
  2750. int diffusion()
  2751. {
  2752.    int xmax,ymax,xmin,ymin;     /* Current maximum coordinates */
  2753.    int border;   /* Distance between release point and fractal */
  2754.    int mode;     /* Determines diffusion type:  0 = central (classic) */
  2755.                  /*                             1 = falling particles */
  2756.                  /*                             2 = square cavity     */
  2757.    int i;
  2758.    double cosine,sine,angle;
  2759.    long lcosine,lsine;
  2760.    int x,y;
  2761.    float r, radius, f_tmp;
  2762.    extern char floatflag;
  2763.  
  2764.    if (diskvideo)
  2765.    {
  2766.       notdiskmsg();
  2767.       return(-1);
  2768.    }
  2769.  
  2770.    bitshift = 16;
  2771.    fudge = 1L << 16;
  2772.  
  2773.    border = param[0];
  2774.    mode = param[1];
  2775.    if (mode > 2)
  2776.       mode=0;
  2777.  
  2778.    if (border <= 0)
  2779.       border = 10;
  2780.  
  2781.    srand(rseed);
  2782.    if (!rflag) ++rseed;
  2783.  
  2784.    if (mode == 0) {
  2785.       xmax = xdots / 2 + border;  /* Initial box */
  2786.       xmin = xdots / 2 - border;
  2787.       ymax = ydots / 2 + border;
  2788.       ymin = ydots / 2 - border;
  2789.    }
  2790.    if (mode == 1) {
  2791.       xmax = xdots / 2 + border;  /* Initial box */
  2792.       xmin = xdots / 2 - border;
  2793.       ymin = ydots - 20;
  2794.    }
  2795.    if (mode == 2) {
  2796.       if (xdots>ydots)
  2797.          radius = ydots - 5;
  2798.       else
  2799.          radius = xdots - 5;
  2800.    }
  2801.    if (resuming) /* restore worklist, if we can't the above will stay in place */
  2802.    {
  2803.       start_resume();
  2804.       if (mode != 2)
  2805.          get_resume(sizeof(int),&xmax,sizeof(int),&xmin,sizeof(int),&ymax,
  2806.              sizeof(int),&ymin,0);
  2807.       else
  2808.          get_resume(sizeof(int),&xmax,sizeof(int),&xmin,sizeof(int),&ymax,
  2809.              sizeof(int),&radius,0);
  2810.       end_resume();
  2811.    }
  2812.  
  2813.    if (mode==0)
  2814.       putcolor(xdots / 2, ydots / 2,RANDOM(colors-1)+1);  /* Seed point */
  2815.  
  2816.    if (mode==1)
  2817.       for (i=0;i<=xdots;i++)
  2818.          putcolor(i,ydots-1,colors-1);
  2819.  
  2820.    if (mode==2){
  2821.       if (xdots>ydots){
  2822.          for (i=0;i<ydots;i++){
  2823.             putcolor(xdots/2-ydots/2 , i , colors-1);
  2824.             putcolor(xdots/2+ydots/2 , i , colors-1);
  2825.             putcolor(xdots/2-ydots/2+i , 0 , colors-1);
  2826.             putcolor(xdots/2-ydots/2+i , ydots-1 , colors-1);
  2827.          }
  2828.       }
  2829.       else {
  2830.          for (i=0;i<xdots;i++){
  2831.             putcolor(0 , ydots/2-xdots/2+i , colors-1);
  2832.             putcolor(xdots-1 , ydots/2-xdots/2+i , colors-1);
  2833.             putcolor(i , ydots/2-xdots/2 , colors-1);
  2834.             putcolor(i , ydots/2+xdots/2 , colors-1);
  2835.          }
  2836.       }
  2837.    }
  2838.  
  2839.    while (1)
  2840.    {
  2841.       /* Release new point on circle just inside the box */
  2842.  
  2843.       if (mode==0)
  2844.       {
  2845.          if (floatflag)
  2846.          {
  2847.             angle=2*(double)rand()/(RAND_MAX/PI);
  2848.             FPUsincos(&angle,&sine,&cosine);
  2849.             x = cosine*(xmax-xmin) + xdots;
  2850.             y = sine  *(ymax-ymin) + ydots;
  2851.          }
  2852.          else
  2853.          {
  2854.             SinCos086(multiply((long)rand15(),FOURPI,16),&lsine,&lcosine);
  2855.             x = (lcosine*(long)(xmax-xmin) >> 16) + xdots;
  2856.             y = (lsine  *(long)(ymax-ymin) >> 16) + ydots;
  2857.          }
  2858.  
  2859.          x = x >> 1; /* divide by 2 */
  2860.          y = y >> 1;
  2861.       }
  2862.       if (mode==1) 
  2863.       {
  2864.          y=ymin;
  2865.          x=RANDOM(xmax-xmin) + (xdots-xmax+xmin)/2;
  2866.       }
  2867.       if (mode==2)
  2868.       {
  2869.          if (floatflag)
  2870.          {
  2871.             angle=2*(double)rand()/(RAND_MAX/PI);
  2872.             FPUsincos(&angle,&sine,&cosine);
  2873.             x = cosine*radius + xdots;
  2874.             y = sine  *radius + ydots;
  2875.          }
  2876.          else
  2877.          {
  2878.             SinCos086(multiply((long)rand15(),FOURPI,16),&lsine,&lcosine);
  2879.             x = (lcosine*(long)(radius) >> 16) + xdots;
  2880.             y = (lsine  *(long)(radius) >> 16) + ydots;
  2881.          }
  2882.          x = x >> 1;
  2883.          y = y >> 1;
  2884.       }
  2885.       /* Loop as long as the point (x,y) is surrounded by color 0 */
  2886.       /* on all eight sides                                       */
  2887.       while((getcolor(x+1,y+1) == 0) && (getcolor(x+1,y) == 0) &&
  2888.           (getcolor(x+1,y-1) == 0) && (getcolor(x  ,y+1) == 0) &&
  2889.           (getcolor(x  ,y-1) == 0) && (getcolor(x-1,y+1) == 0) &&
  2890.           (getcolor(x-1,y) == 0) && (getcolor(x-1,y-1) == 0))
  2891.       {
  2892.          /* Erase moving point */
  2893.          if (show_orbit)
  2894.             putcolor(x,y,0);
  2895.  
  2896.          /* Make sure point is inside the box (if mode==0)*/
  2897.          if (mode==0){
  2898.             if (x==xmax)
  2899.                x--;
  2900.             else if (x==xmin)
  2901.                x++;
  2902.             if (y==ymax)
  2903.                y--;
  2904.             else if (y==ymin)
  2905.                y++;
  2906.          }
  2907.  
  2908.          if (mode==1) 
  2909.          {
  2910.             if (x>xdots-2)
  2911.                x--;
  2912.             else if (x<1)
  2913.                x++;
  2914.             if (y<ymin)
  2915.                y++;
  2916.          }
  2917.  
  2918.          /* Take one random step */
  2919.          x += RANDOM(3) - 1;
  2920.          y += RANDOM(3) - 1;
  2921.  
  2922.          /* Check keyboard */
  2923.          if ((++plasma_check & 0x7f) == 1)
  2924.             if(check_key())
  2925.             {
  2926.                alloc_resume(20,1);
  2927.                if (mode!=2)
  2928.                   put_resume(sizeof(int),&xmax,sizeof(int),&xmin, sizeof(int),&ymax,
  2929.                       sizeof(int),&ymin,0);
  2930.                else
  2931.                   put_resume(sizeof(int),&xmax,sizeof(int),&xmin, sizeof(int),&ymax,
  2932.                       sizeof(int),&radius,0);
  2933.  
  2934.                plasma_check--;
  2935.                return 1;
  2936.             }
  2937.  
  2938.          /* Show the moving point */
  2939.          if (show_orbit)
  2940.             putcolor(x,y,RANDOM(colors-1)+1);
  2941.  
  2942.       }
  2943.       putcolor(x,y,RANDOM(colors-1)+1);
  2944.  
  2945.       /* Is point too close to the edge? */
  2946.  
  2947.       if (mode==0)
  2948.       {
  2949.          if (((x+border)>xmax) || ((x-border)<xmin)
  2950.              || ((y-border)<ymin) || ((y+border)>ymax))
  2951.          {
  2952.             /* Increase box size, but not past the edge of the screen */
  2953.             if (ymin != 1)
  2954.             {
  2955.                ymin--;
  2956.                ymax++;
  2957.             }
  2958.             if (xmin != 1)
  2959.             {
  2960.                xmin--;
  2961.                xmax++;
  2962.             }
  2963.             if ((ymin==1) || (xmin==1))
  2964.                return 0;
  2965.          }
  2966.       }
  2967.       if (mode==1) 
  2968.       {
  2969.          if (y < ymin+5)
  2970.             ymin = y - 5;
  2971.          if (ymin<2)
  2972.             return 0;
  2973.       }
  2974.       if (mode==2)
  2975.       {
  2976.          if (abs(x-xdots/2)<5 && abs(y-ydots/2)<5)
  2977.             return 0;
  2978.  
  2979.          r = (x-xdots/2)*(x-xdots/2)+(y-ydots/2)*(y-ydots/2);
  2980.          fSqrt14(r,f_tmp);
  2981.          r = 2 * f_tmp;
  2982.          if (r < radius)
  2983.             radius = r;
  2984.       }
  2985.    }
  2986. }
  2987.  
  2988. /************ standalone engine for "bifurcation" types ***************/
  2989.  
  2990. /***************************************************************/
  2991. /* The following code now forms a generalised Fractal Engine   */
  2992. /* for Bifurcation fractal typeS.  By rights it now belongs in */
  2993. /* CALCFRACT.C, but it's easier for me to leave it here !      */
  2994.  
  2995. /* Original code by Phil Wilson, hacked around by Kev Allen.   */
  2996.  
  2997. /* Besides generalisation, enhancements include Periodicity    */
  2998. /* Checking during the plotting phase (AND halfway through the */
  2999. /* filter cycle, if possible, to halve calc times), quicker    */
  3000. /* floating-point calculations for the standard Verhulst type, */
  3001. /* and new bifurcation types (integer bifurcation, f.p & int   */
  3002. /* biflambda - the real equivalent of complex Lambda sets -    */
  3003. /* and f.p renditions of bifurcations of r*sin(Pi*p), which    */
  3004. /* spurred Mitchel Feigenbaum on to discover his Number).      */
  3005.  
  3006. /* To add further types, extend the fractalspecific[] array in */
  3007. /* usual way, with Bifurcation as the engine, and the name of  */
  3008. /* the routine that calculates the next bifurcation generation */
  3009. /* as the "orbitcalc" routine in the fractalspecific[] entry.  */
  3010.  
  3011. /* Bifurcation "orbitcalc" routines get called once per screen */
  3012. /* pixel column.  They should calculate the next generation    */
  3013. /* from the doubles Rate & Population (or the longs lRate &    */
  3014. /* lPopulation if they use integer math), placing the result   */
  3015. /* back in Population (or lPopulation).  They should return 0  */
  3016. /* if all is ok, or any non-zero value if calculation bailout  */
  3017. /* is desirable (eg in case of errors, or the series tending   */
  3018. /* to infinity).        Have fun !               */
  3019. /***************************************************************/
  3020.  
  3021. #define BIG 100000.0
  3022. #define DEFAULTFILTER 1000     /* "Beauty of Fractals" recommends using 5000
  3023.                    (p.25), but that seems unnecessary. Can
  3024.                    override this value with a nonzero param1 */
  3025.  
  3026. #define SEED 0.66        /* starting value for population */
  3027.  
  3028. static int far *verhulst_array;
  3029. static unsigned int filter_cycles;
  3030. static unsigned int half_time_check;
  3031. static long   lPopulation, lRate;
  3032. double Population,  Rate;
  3033. static int    mono, outside_x;
  3034. static long   LPI;
  3035.  
  3036. int Bifurcation(void)
  3037. {
  3038.    unsigned long array_size;
  3039.    int row, column;
  3040.    column = 0;
  3041.    if (resuming)
  3042.    {
  3043.       start_resume();
  3044.       get_resume(sizeof(int),&column,0);
  3045.       end_resume();
  3046.    }
  3047.    array_size =  (iystop + 1) * sizeof(int); /* should be iystop + 1 */
  3048.    if ((verhulst_array = (int far *) farmemalloc(array_size)) == NULL)
  3049.    {
  3050.       static char far msg[]={"Insufficient free memory for calculation."};
  3051.       stopmsg(0,msg);
  3052.       return(-1);
  3053.    }
  3054.  
  3055.    LPI = PI * fudge;
  3056.  
  3057.    for (row = 0; row <= iystop; row++) /* should be iystop */
  3058.       verhulst_array[row] = 0;
  3059.  
  3060.    mono = 0;
  3061.    if (colors == 2)
  3062.       mono = 1;
  3063.    if (mono)
  3064.    {
  3065.       if (inside)
  3066.       {
  3067.      outside_x = 0;
  3068.      inside = 1;
  3069.       }
  3070.       else
  3071.      outside_x = 1;
  3072.    }
  3073.  
  3074.    filter_cycles = (parm.x <= 0) ? DEFAULTFILTER : parm.x;
  3075.    half_time_check = FALSE;
  3076.    if (periodicitycheck && maxit < filter_cycles)
  3077.    {
  3078.       filter_cycles = (filter_cycles - maxit + 1) / 2;
  3079.       half_time_check = TRUE;
  3080.    }
  3081.  
  3082.    if (integerfractal)
  3083.       linit.y = ly0[iystop];   /* Y-value of    */
  3084.    else
  3085.       init.y = dy0[iystop];   /* bottom pixels */
  3086.  
  3087.    while (column <= ixstop)
  3088.    {
  3089.       if(check_key())
  3090.       {
  3091.      farmemfree((char far *)verhulst_array);
  3092.      alloc_resume(10,1);
  3093.      put_resume(sizeof(int),&column,0);
  3094.      return(-1);
  3095.       }
  3096.       if (integerfractal)
  3097.      lRate = lx0[column];
  3098.       else
  3099.      Rate = dx0[column];
  3100.       verhulst();     /* calculate array once per column */
  3101.  
  3102.       for (row = iystop; row >= 0; row--) /* should be iystop & >=0 */
  3103.       {
  3104.      int color;
  3105.      color = verhulst_array[row];
  3106.      if(color && mono)
  3107.         color = inside;
  3108.      else if((!color) && mono)
  3109.         color = outside_x;
  3110.      else if (color>=colors)
  3111.         color = colors-1;
  3112.      verhulst_array[row] = 0;
  3113.      (*plot)(column,row,color); /* was row-1, but that's not right? */
  3114.       }
  3115.       column++;
  3116.    }
  3117.    farmemfree((char far *)verhulst_array);
  3118.    return(0);
  3119. }
  3120.  
  3121. static void verhulst()        /* P. F. Verhulst (1845) */
  3122. {
  3123.    unsigned int pixel_row, counter, errors;
  3124.  
  3125.     if (integerfractal)
  3126.        lPopulation = (parm.y == 0) ? SEED * fudge : parm.y * fudge;
  3127.     else
  3128.        Population = (parm.y == 0 ) ? SEED : parm.y;
  3129.  
  3130.    errors = overflow = FALSE;
  3131.  
  3132.    for (counter=0 ; counter < filter_cycles ; counter++)
  3133.    {
  3134.       errors = (*(curfractalspecific->orbitcalc))();
  3135.       if (errors)
  3136.      return;
  3137.    }
  3138.    if (half_time_check) /* check for periodicity at half-time */
  3139.    {
  3140.       Bif_Period_Init();
  3141.       for (counter=0 ; counter < maxit ; counter++)
  3142.       {
  3143.      errors = (*(curfractalspecific->orbitcalc))();
  3144.      if (errors) return;
  3145.      if (periodicitycheck && Bif_Periodic(counter)) break;
  3146.       }
  3147.       if (counter >= maxit)   /* if not periodic, go the distance */
  3148.       {
  3149.      for (counter=0 ; counter < filter_cycles ; counter++)
  3150.      {
  3151.         errors = (*(curfractalspecific->orbitcalc))();
  3152.         if (errors) return;
  3153.      }
  3154.       }
  3155.    }
  3156.  
  3157.    if (periodicitycheck) Bif_Period_Init();
  3158.    for (counter=0 ; counter < maxit ; counter++)
  3159.    {
  3160.       errors = (*(curfractalspecific->orbitcalc))();
  3161.       if (errors) return;
  3162.  
  3163.       /* assign population value to Y coordinate in pixels */
  3164.       if (integerfractal)
  3165.      pixel_row = iystop - (lPopulation - linit.y) / dely; /* iystop */
  3166.       else
  3167.      pixel_row = iystop - (int)((Population - init.y) / deltaY);
  3168.  
  3169.       /* if it's visible on the screen, save it in the column array */
  3170.       if (pixel_row <= iystop) /* JCO 6/6/92 */
  3171.      verhulst_array[ pixel_row ] ++;
  3172.       if (periodicitycheck && Bif_Periodic(counter))
  3173.       {
  3174.      if (pixel_row <= iystop) /* JCO 6/6/92 */
  3175.         verhulst_array[ pixel_row ] --;
  3176.      break;
  3177.       }
  3178.    }
  3179. }
  3180. static    long    lBif_closenuf, lBif_savedpop;    /* poss future use  */
  3181. static    double     Bif_closenuf,    Bif_savedpop;
  3182. static    int     Bif_savedinc,    Bif_savedand;
  3183.  
  3184. static void Bif_Period_Init()
  3185. {
  3186.    Bif_savedinc = 1;
  3187.    Bif_savedand = 1;
  3188.    if (integerfractal)
  3189.    {
  3190.       lBif_savedpop = -1;
  3191.       lBif_closenuf = dely / 8;
  3192.    }
  3193.    else
  3194.    {
  3195.       Bif_savedpop = -1.0;
  3196.       Bif_closenuf = deltaY / 8.0;
  3197.    }
  3198. }
  3199.  
  3200. static int _fastcall Bif_Periodic (time)  /* Bifurcation Population Periodicity Check */
  3201. int time;        /* Returns : 1 if periodicity found, else 0 */
  3202. {
  3203.    if ((time & Bif_savedand) == 0)    /* time to save a new value */
  3204.    {
  3205.       if (integerfractal) lBif_savedpop = lPopulation;
  3206.       else             Bif_savedpop =  Population;
  3207.       if (--Bif_savedinc == 0)
  3208.       {
  3209.      Bif_savedand = (Bif_savedand << 1) + 1;
  3210.      Bif_savedinc = 4;
  3211.       }
  3212.    }
  3213.    else             /* check against an old save */
  3214.    {
  3215.       if (integerfractal)
  3216.       {
  3217.      if (labs(lBif_savedpop-lPopulation) <= lBif_closenuf)
  3218.         return(1);
  3219.       }
  3220.       else
  3221.       {
  3222.      if (fabs(Bif_savedpop-Population) <= Bif_closenuf)
  3223.         return(1);
  3224.       }
  3225.    }
  3226.    return(0);
  3227. }
  3228.  
  3229. /**********************************************************************/
  3230. /*                                                      */
  3231. /* The following are Bifurcation "orbitcalc" routines...              */
  3232. /*                                                      */
  3233. /**********************************************************************/
  3234. #ifdef XFRACT
  3235. int BifurcLambda() /* Used by lyanupov */
  3236.   {
  3237.     Population = Rate * Population * (1 - Population);
  3238.     return (fabs(Population) > BIG);
  3239.   }
  3240. #endif
  3241.  
  3242. /* Modified formulas below to generalize bifurcations. JCO 7/3/92 */
  3243.  
  3244. #define LCMPLXtrig0(arg,out) Arg1->l = (arg); ltrig0(); (out)=Arg1->l
  3245. #define  CMPLXtrig0(arg,out) Arg1->d = (arg); dtrig0(); (out)=Arg1->d
  3246.  
  3247. int BifurcVerhulstTrig()
  3248.   {
  3249. /*  Population = Pop + Rate * fn(Pop) * (1 - fn(Pop)) */
  3250.     tmp.x = Population;
  3251.     tmp.y = 0;
  3252.     CMPLXtrig0(tmp, tmp);
  3253.     Population += Rate * tmp.x * (1 - tmp.x);
  3254.     return (fabs(Population) > BIG);
  3255.   }
  3256.  
  3257. int LongBifurcVerhulstTrig()
  3258.   {
  3259. #ifndef XFRACT
  3260.     ltmp.x = lPopulation;
  3261.     ltmp.y = 0;
  3262.     LCMPLXtrig0(ltmp, ltmp);
  3263.     ltmp.y = ltmp.x - multiply(ltmp.x,ltmp.x,bitshift);
  3264.     lPopulation += multiply(lRate,ltmp.y,bitshift);
  3265.     return (overflow);
  3266. #endif
  3267.   }
  3268.  
  3269. int BifurcStewartTrig()
  3270.   {
  3271. /*  Population = (Rate * fn(Population) * fn(Population)) - 1.0 */
  3272.     tmp.x = Population;
  3273.     tmp.y = 0;
  3274.     CMPLXtrig0(tmp, tmp);
  3275.     Population = (Rate * tmp.x * tmp.x) - 1.0;
  3276.     return (fabs(Population) > BIG);
  3277.   }
  3278.  
  3279. int LongBifurcStewartTrig()
  3280.   {
  3281. #ifndef XFRACT
  3282.     ltmp.x = lPopulation;
  3283.     ltmp.y = 0;
  3284.     LCMPLXtrig0(ltmp, ltmp);
  3285.     lPopulation = multiply(ltmp.x,ltmp.x,bitshift);
  3286.     lPopulation = multiply(lPopulation,lRate,       bitshift);
  3287.     lPopulation -= fudge;
  3288.     return (overflow);
  3289. #endif
  3290.   }
  3291.  
  3292. int BifurcSetTrigPi()
  3293.   {
  3294.     tmp.x = Population * PI;
  3295.     tmp.y = 0;
  3296.     CMPLXtrig0(tmp, tmp);
  3297.     Population = Rate * tmp.x;
  3298.     return (fabs(Population) > BIG);
  3299.   }
  3300.  
  3301. int LongBifurcSetTrigPi()
  3302.   {
  3303. #ifndef XFRACT
  3304.     ltmp.x = multiply(lPopulation,LPI,bitshift);
  3305.     ltmp.y = 0;
  3306.     LCMPLXtrig0(ltmp, ltmp);
  3307.     lPopulation = multiply(lRate,ltmp.x,bitshift);
  3308.     return (overflow);
  3309. #endif
  3310.   }
  3311.  
  3312. int BifurcAddTrigPi()
  3313.   {
  3314.     tmp.x = Population * PI;
  3315.     tmp.y = 0;
  3316.     CMPLXtrig0(tmp, tmp);
  3317.     Population += Rate * tmp.x;
  3318.     return (fabs(Population) > BIG);
  3319.   }
  3320.  
  3321. int LongBifurcAddTrigPi()
  3322.   {
  3323. #ifndef XFRACT
  3324.     ltmp.x = multiply(lPopulation,LPI,bitshift);
  3325.     ltmp.y = 0;
  3326.     LCMPLXtrig0(ltmp, ltmp);
  3327.     lPopulation += multiply(lRate,ltmp.x,bitshift);
  3328.     return (overflow);
  3329. #endif
  3330.   }
  3331.  
  3332. int BifurcLambdaTrig()
  3333.   {
  3334. /*  Population = Rate * fn(Population) * (1 - fn(Population)) */
  3335.     tmp.x = Population;
  3336.     tmp.y = 0;
  3337.     CMPLXtrig0(tmp, tmp);
  3338.     Population = Rate * tmp.x * (1 - tmp.x);
  3339.     return (fabs(Population) > BIG);
  3340.   }
  3341.  
  3342. int LongBifurcLambdaTrig()
  3343.   {
  3344. #ifndef XFRACT
  3345.     ltmp.x = lPopulation;
  3346.     ltmp.y = 0;
  3347.     LCMPLXtrig0(ltmp, ltmp);
  3348.     ltmp.y = ltmp.x - multiply(ltmp.x,ltmp.x,bitshift);
  3349.     lPopulation = multiply(lRate,ltmp.y,bitshift);
  3350.     return (overflow);
  3351. #endif
  3352.   }
  3353.  
  3354. #define LCMPLXpwr(arg1,arg2,out)    Arg2->l = (arg1); Arg1->l = (arg2);\
  3355.      lStkPwr(); Arg1++; Arg2++; (out) = Arg2->l
  3356.  
  3357. int BifurcMay()
  3358.   { /* X = (lambda * X) / (1 + X)^beta, from R.May as described in Pickover,
  3359.             Computers, Pattern, Chaos, and Beauty, page 153 */
  3360.     tmp.x = 1.0 + Population;
  3361.     tmp.x = pow(tmp.x, -parm2.x); /* pow in math.h included with mpmath.h */
  3362.     Population = (Rate * Population) * tmp.x;
  3363.     return (fabs(Population) > BIG);
  3364.   }
  3365.  
  3366. int LongBifurcMay()
  3367.   {
  3368. #ifndef XFRACT
  3369.     ltmp.x = lPopulation + fudge;
  3370.     ltmp.y = 0;
  3371.     LCMPLXpwr(ltmp, lparm2, ltmp);
  3372.     lPopulation = multiply(lRate,lPopulation,bitshift);
  3373.     lPopulation = divide(lPopulation,ltmp.x,bitshift);
  3374.     return (overflow);
  3375. #endif
  3376.   }
  3377.  
  3378.  
  3379. /* Here Endeth the Generalised Bifurcation Fractal Engine   */
  3380.  
  3381. /* END Phil Wilson's Code (modified slightly by Kev Allen et. al. !) */
  3382.  
  3383.  
  3384. /******************* standalone engine for "popcorn" ********************/
  3385.  
  3386. int popcorn()    /* subset of std engine */
  3387. {
  3388.    int start_row;
  3389.    start_row = 0;
  3390.    if (resuming)
  3391.    {
  3392.       start_resume();
  3393.       get_resume(sizeof(int),&start_row,0);
  3394.       end_resume();
  3395.    }
  3396.    kbdcount=max_kbdcount;
  3397.    plot = noplot;
  3398.    tempsqrx = ltempsqrx = 0; /* PB added this to cover weird BAILOUTs */
  3399.    for (row = start_row; row <= iystop; row++)
  3400.    {
  3401.       reset_periodicity = 1;
  3402.       for (col = 0; col <= ixstop; col++)
  3403.       {
  3404.      if (StandardFractal() == -1) /* interrupted */
  3405.      {
  3406.         alloc_resume(10,1);
  3407.         put_resume(sizeof(int),&row,0);
  3408.         return(-1);
  3409.      }
  3410.      reset_periodicity = 0;
  3411.       }
  3412.    }
  3413.    calc_status = 4;
  3414.    return(0);
  3415. }
  3416.  
  3417. /**************** tesseral method by CJLT begins here*********************/
  3418. /*  reworked by PB for speed and resumeability */
  3419.  
  3420. struct tess { /* one of these per box to be done gets stacked */
  3421.    int x1,x2,y1,y2;     /* left/right top/bottom x/y coords  */
  3422.    int top,bot,lft,rgt;  /* edge colors, -1 mixed, -2 unknown */
  3423. };
  3424.  
  3425. static int tesseral()
  3426. {
  3427.    register struct tess *tp;
  3428.  
  3429.    guessplot = (plot != putcolor && plot != symplot2);
  3430.    tp = (struct tess *)&dstack[0];
  3431.    tp->x1 = ixstart;                  /* set up initial box */
  3432.    tp->x2 = ixstop;
  3433.    tp->y1 = iystart;
  3434.    tp->y2 = iystop;
  3435.  
  3436.    if (workpass == 0) { /* not resuming */
  3437.       tp->top = tessrow(ixstart,ixstop,iystart);     /* Do top row */
  3438.       tp->bot = tessrow(ixstart,ixstop,iystop);      /* Do bottom row */
  3439.       tp->lft = tesscol(ixstart,iystart+1,iystop-1); /* Do left column */
  3440.       tp->rgt = tesscol(ixstop,iystart+1,iystop-1);  /* Do right column */
  3441.       if (check_key()) { /* interrupt before we got properly rolling */
  3442.          add_worklist(xxstart,xxstop,yystart,yystop,yystart,0,worksym);
  3443.          return(-1);
  3444.       }
  3445.    }
  3446.  
  3447.    else { /* resuming, rebuild work stack */
  3448.       int i,mid,curx,cury,xsize,ysize;
  3449.       struct tess *tp2;
  3450.       tp->top = tp->bot = tp->lft = tp->rgt = -2;
  3451.       cury = yybegin & 0xfff;
  3452.       ysize = 1;
  3453.       i = (unsigned)yybegin >> 12;
  3454.       while (--i >= 0) ysize <<= 1;
  3455.       curx = workpass & 0xfff;
  3456.       xsize = 1;
  3457.       i = (unsigned)workpass >> 12;
  3458.       while (--i >= 0) xsize <<= 1;
  3459.       while (1) {
  3460.          tp2 = tp;
  3461.          if (tp->x2 - tp->x1 > tp->y2 - tp->y1) { /* next divide down middle */
  3462.             if (tp->x1 == curx && (tp->x2 - tp->x1 - 2) < xsize)
  3463.                break;
  3464.             mid = (tp->x1 + tp->x2) >> 1;         /* Find mid point */
  3465.             if (mid > curx) { /* stack right part */
  3466.                memcpy(++tp,tp2,sizeof(*tp));
  3467.                tp->x2 = mid;
  3468.             }
  3469.             tp2->x1 = mid;
  3470.          }
  3471.          else {                   /* next divide across */
  3472.             if (tp->y1 == cury && (tp->y2 - tp->y1 - 2) < ysize) break;
  3473.             mid = (tp->y1 + tp->y2) >> 1;         /* Find mid point */
  3474.             if (mid > cury) { /* stack bottom part */
  3475.                memcpy(++tp,tp2,sizeof(*tp));
  3476.                tp->y2 = mid;
  3477.             }
  3478.             tp2->y1 = mid;
  3479.          }
  3480.       }
  3481.    }
  3482.  
  3483.    got_status = 4; /* for tab_display */
  3484.  
  3485.    while (tp >= (struct tess *)&dstack[0]) { /* do next box */
  3486.       curcol = tp->x1; /* for tab_display */
  3487.       currow = tp->y1;
  3488.  
  3489.       if (tp->top == -1 || tp->bot == -1 || tp->lft == -1 || tp->rgt == -1)
  3490.          goto tess_split;
  3491.       /* for any edge whose color is unknown, set it */
  3492.       if (tp->top == -2)
  3493.          tp->top = tesschkrow(tp->x1,tp->x2,tp->y1);
  3494.       if (tp->top == -1)
  3495.          goto tess_split;
  3496.       if (tp->bot == -2)
  3497.          tp->bot = tesschkrow(tp->x1,tp->x2,tp->y2);
  3498.       if (tp->bot != tp->top)
  3499.          goto tess_split;
  3500.       if (tp->lft == -2)
  3501.          tp->lft = tesschkcol(tp->x1,tp->y1,tp->y2);
  3502.       if (tp->lft != tp->top)
  3503.          goto tess_split;
  3504.       if (tp->rgt == -2)
  3505.          tp->rgt = tesschkcol(tp->x2,tp->y1,tp->y2);
  3506.       if (tp->rgt != tp->top)
  3507.          goto tess_split;
  3508.  
  3509.       {  /* all 4 edges are the same color, fill in */
  3510.          int i,j;
  3511.          i = 0;
  3512.          if(fillcolor != 0)
  3513.          {
  3514.          if(fillcolor > 0)
  3515.             tp->top = fillcolor & (colors-1);
  3516.          if (guessplot || (j = tp->x2 - tp->x1 - 1) < 2) { /* paint dots */
  3517.             for (col = tp->x1 + 1; col < tp->x2; col++)
  3518.                for (row = tp->y1 + 1; row < tp->y2; row++) {
  3519.                   (*plot)(col,row,tp->top);
  3520.                   if (++i > 500) {
  3521.                      if (check_key()) goto tess_end;
  3522.                      i = 0;
  3523.                   }
  3524.                }
  3525.          }
  3526.          else { /* use put_line for speed */
  3527.             memset(&dstack[MAXPIXELS],tp->top,j);
  3528.             for (row = tp->y1 + 1; row < tp->y2; row++) {
  3529.                put_line(row,tp->x1+1,tp->x2-1,&dstack[MAXPIXELS]);
  3530.                if (plot != putcolor) /* symmetry */
  3531.                   if ((j = yystop-(row-yystart)) > iystop && j < ydots)
  3532.                      put_line(j,tp->x1+1,tp->x2-1,&dstack[MAXPIXELS]);
  3533.                if (++i > 25) {
  3534.                   if (check_key()) goto tess_end;
  3535.                   i = 0;
  3536.                }
  3537.             }
  3538.          }
  3539.          }
  3540.          --tp;
  3541.       }
  3542.       continue;
  3543.  
  3544.       tess_split:
  3545.       {  /* box not surrounded by same color, sub-divide */
  3546.          int mid,midcolor;
  3547.          struct tess *tp2;
  3548.          if (tp->x2 - tp->x1 > tp->y2 - tp->y1) { /* divide down the middle */
  3549.             mid = (tp->x1 + tp->x2) >> 1;         /* Find mid point */
  3550.             midcolor = tesscol(mid, tp->y1+1, tp->y2-1); /* Do mid column */
  3551.             if (midcolor == -3) goto tess_end;
  3552.             if (tp->x2 - mid > 1) {    /* right part >= 1 column */
  3553.                if (tp->top == -1) tp->top = -2;
  3554.                if (tp->bot == -1) tp->bot = -2;
  3555.                tp2 = tp;
  3556.                if (mid - tp->x1 > 1) { /* left part >= 1 col, stack right */
  3557.                   memcpy(++tp,tp2,sizeof(*tp));
  3558.                   tp->x2 = mid;
  3559.                   tp->rgt = midcolor;
  3560.                }
  3561.                tp2->x1 = mid;
  3562.                tp2->lft = midcolor;
  3563.             }
  3564.             else
  3565.                --tp;
  3566.          }
  3567.          else {                   /* divide across the middle */
  3568.             mid = (tp->y1 + tp->y2) >> 1;         /* Find mid point */
  3569.             midcolor = tessrow(tp->x1+1, tp->x2-1, mid); /* Do mid row */
  3570.             if (midcolor == -3) goto tess_end;
  3571.             if (tp->y2 - mid > 1) {    /* bottom part >= 1 column */
  3572.                if (tp->lft == -1) tp->lft = -2;
  3573.                if (tp->rgt == -1) tp->rgt = -2;
  3574.                tp2 = tp;
  3575.                if (mid - tp->y1 > 1) { /* top also >= 1 col, stack bottom */
  3576.                   memcpy(++tp,tp2,sizeof(*tp));
  3577.                   tp->y2 = mid;
  3578.                   tp->bot = midcolor;
  3579.                }
  3580.                tp2->y1 = mid;
  3581.                tp2->top = midcolor;
  3582.             }
  3583.             else
  3584.                --tp;
  3585.          }
  3586.       }
  3587.  
  3588.    }
  3589.  
  3590.    tess_end:
  3591.    if (tp >= (struct tess *)&dstack[0]) { /* didn't complete */
  3592.       int i,xsize,ysize;
  3593.       xsize = ysize = 1;
  3594.       i = 2;
  3595.       while (tp->x2 - tp->x1 - 2 >= i) {
  3596.          i <<= 1;
  3597.          ++xsize;
  3598.       }
  3599.       i = 2;
  3600.       while (tp->y2 - tp->y1 - 2 >= i) {
  3601.          i <<= 1;
  3602.          ++ysize;
  3603.       }
  3604.       add_worklist(xxstart,xxstop,yystart,yystop,
  3605.           (ysize<<12)+tp->y1,(xsize<<12)+tp->x1,worksym);
  3606.       return(-1);
  3607.    }
  3608.    return(0);
  3609.  
  3610. } /* tesseral */
  3611.  
  3612. static int _fastcall tesschkcol(int x,int y1,int y2)
  3613. {
  3614.    int i;
  3615.    i = getcolor(x,++y1);
  3616.    while (--y2 > y1)
  3617.       if (getcolor(x,y2) != i) return -1;
  3618.    return i;
  3619. }
  3620.  
  3621. static int _fastcall tesschkrow(int x1,int x2,int y)
  3622. {
  3623.    int i;
  3624.    i = getcolor(x1,y);
  3625.    while (x2 > x1) {
  3626.       if (getcolor(x2,y) != i) return -1;
  3627.       --x2;
  3628.    }
  3629.    return i;
  3630. }
  3631.  
  3632. static int _fastcall tesscol(int x,int y1,int y2)
  3633. {
  3634.    int colcolor,i;
  3635.    col = x;
  3636.    row = y1;
  3637.    reset_periodicity = 1;
  3638.    colcolor = (*calctype)();
  3639.    reset_periodicity = 0;
  3640.    while (++row <= y2) { /* generate the column */
  3641.       if ((i = (*calctype)()) < 0) return -3;
  3642.       if (i != colcolor) colcolor = -1;
  3643.    }
  3644.    return colcolor;
  3645. }
  3646.  
  3647. static int _fastcall tessrow(int x1,int x2,int y)
  3648. {
  3649.    int rowcolor,i;
  3650.    row = y;
  3651.    col = x1;
  3652.    reset_periodicity = 1;
  3653.    rowcolor = (*calctype)();
  3654.    reset_periodicity = 0;
  3655.    while (++col <= x2) { /* generate the row */
  3656.       if ((i = (*calctype)()) < 0) return -3;
  3657.       if (i != rowcolor) rowcolor = -1;
  3658.    }
  3659.    return rowcolor;
  3660. }
  3661.  
  3662. /******************* standalone engine for "lyapunov" ********************/
  3663. /*** Roy Murphy [76376,721] ***/
  3664. /*** benchmark 15:05        ***/
  3665. int lyaLength;
  3666. int lyaRxy[18];
  3667.  
  3668. extern int lyapunov_cycles(int, double, double, double);
  3669.  
  3670. lyapunov () {
  3671.     int count, i, row, col, startrow, lnadjust, filter_cycles;
  3672.     double a, b, lyap, total, temp;
  3673.     /* e10=22026.4657948  e-10=0.0000453999297625 */
  3674.  
  3675.       startrow = 0;
  3676.       if (resuming) {
  3677.           start_resume();
  3678.           get_resume(sizeof(int), &startrow, 0);
  3679.           end_resume();
  3680.           }
  3681.       for (row=startrow; row<iystop; row++) {
  3682.            for (col=0; col<ixstop; col++) {
  3683.             if (check_key()) {
  3684.                 alloc_resume(10,1);
  3685.                 put_resume(sizeof(int), &row, 0);
  3686.         return -1;
  3687.                 }
  3688.  
  3689.             overflow=FALSE;
  3690.             total = 1.0;
  3691.             lnadjust = 0;
  3692.                (*plot)(col, row, 1);
  3693.  
  3694.             a = dx0[col]+dx1[row];
  3695.             b = dy0[row]+dy1[col];
  3696.             Population = param[1];
  3697.             if ((filter_cycles=param[2])==0)
  3698.                 filter_cycles=maxit/2;
  3699. #ifndef XFRACT
  3700.         (*plot)(col, row, lyapunov_cycles(filter_cycles,
  3701.                           Population,
  3702.                           a, b));
  3703. #else
  3704.             for (i=0; i<filter_cycles; i++) {
  3705.                 for (count=0; count<lyaLength; count++) {
  3706.             Rate = lyaRxy[count] ? a : b;
  3707.                     if (curfractalspecific->orbitcalc()) {
  3708.                         overflow = TRUE;
  3709.                         goto jumpout;
  3710.                         }
  3711.                     }
  3712.                 }
  3713.  
  3714.         for (i=0; i < maxit/2; i++) {
  3715.         for (count = 0; count < lyaLength; count++) {
  3716.             Rate = lyaRxy[count] ? a : b;
  3717.                     if (curfractalspecific->orbitcalc()) {
  3718.                         overflow = TRUE;
  3719.                         goto jumpout;
  3720.                         }
  3721.                     temp = fabs(Rate-2.0*Rate*Population);
  3722.                     if ((total *= (temp))==0) {
  3723.                         lnadjust = 0;
  3724.                         goto jumpout;
  3725.                         }
  3726.                     while (total > 22026.4657948) {
  3727.                         total *= 0.0000453999297625;
  3728.                         lnadjust += 10;
  3729.                         }
  3730.                     while (total < 0.0000453999297625) {
  3731.                         total *= 22026.4657948;
  3732.                         lnadjust -= 10;
  3733.                         }
  3734.                     }
  3735.                    }
  3736.  
  3737. jumpout:    if (overflow ||
  3738.         total <= 0 ||
  3739.         (temp = log(total) + lnadjust) > 0)
  3740.                 color = 0;
  3741.             else {
  3742.         if (LogFlag)
  3743.                     lyap = -temp/((double) lyaLength*i);
  3744.                 else
  3745.                     lyap = 1 - exp(temp/((double) lyaLength*i));
  3746.                 color = 1 + (int)(lyap * (colors-1));
  3747.                 }
  3748.             (*plot)(col, row, color);
  3749. #endif
  3750.             }
  3751.         }
  3752.     return 1;
  3753. }
  3754.  
  3755. lya_setup () {
  3756.       /* This routine sets up the sequence for forcing the Rate parameter
  3757.           to vary between the two values.  It fills the array lyaRxy[] and
  3758.           sets lyaLength to the length of the sequence.
  3759.  
  3760.           The sequence is coded in the bit pattern in an integer.
  3761.           Briefly, the sequence starts with an A the leading zero bits
  3762.           are ignored and the remaining bit sequence is decoded.  The
  3763.           sequence ends with a B.  Not all possible sequences can be
  3764.           represented in this manner, but every possible sequence is
  3765.           either represented as itself, as a rotation of one of the
  3766.           representable sequences, or as the inverse of a representable
  3767.           sequence (swapping 0s and 1s in the array.)  Sequences that
  3768.           are the rotation and/or inverses of another sequence will generate
  3769.           the same lyapunov exponents.
  3770.  
  3771.           A few examples follow:
  3772.               number    sequence
  3773.                   0        ab
  3774.                   1        aab
  3775.                   2        aabb
  3776.                   3        aaab
  3777.                   4        aabbb
  3778.                   5        aabab
  3779.                   6        aaabb (this is a duplicate of 4, a rotated inverse)
  3780.                   7        aaaab
  3781.                   8        aabbbb    etc.
  3782.      */
  3783.  
  3784.     int i, t;
  3785.  
  3786.  
  3787.     lyaLength = 1;
  3788.     i = (int) param[0];
  3789.     lyaRxy[0] = 1;
  3790.     for (t=15; t>=0; t--)
  3791.         if (i & (1<<t)) break;
  3792.     for (; t>=0; t--)
  3793.         lyaRxy[lyaLength++] = (i & (1<<t)) != 0;
  3794.     lyaRxy[lyaLength++] = 0;
  3795.     return 1;
  3796.  }
  3797.  
  3798. /******************* standalone engine for "cellular" ********************/
  3799. /* Originally coded by Ken Shirriff.
  3800.    Modified beyond recognition by Jonathan Osuch.
  3801.      Original or'd the neighborhood, changed to sum the neighborhood
  3802.      Changed prompts and error messages
  3803.      Added CA types
  3804.      Set the palette to some standard? CA colors
  3805.      Changed *cell_array to near and used dstack so put_line and get_line
  3806.        could be used all the time
  3807.      Made space bar generate next screen
  3808. */
  3809.  
  3810. #define BAD_T         1
  3811. #define BAD_MEM       2
  3812. #define STRING1       3
  3813. #define STRING2       4
  3814. #define TABLEK        5
  3815. #define TYPEKR        6
  3816. #define RULELENGTH    7
  3817.  
  3818. #define CELLULAR_DONE 10
  3819.  
  3820. #ifndef XFRACT
  3821. static BYTE *cell_array[2];
  3822. #else
  3823. static BYTE far *cell_array[2];
  3824. #endif
  3825.  
  3826. S16 r, k_1, rule_digits;
  3827.  
  3828. void abort_cellular(int err, int t)
  3829. {
  3830.    int i;
  3831.    switch (err)
  3832.    {
  3833.       case BAD_T: 
  3834.          {
  3835.          char msg[30];
  3836.          sprintf(msg,"Bad t=%d, aborting\n", t);
  3837.          stopmsg(0,(char far *)msg);
  3838.          }
  3839.          break;
  3840.       case BAD_MEM:
  3841.          {
  3842.          static char far msg[]={"Insufficient free memory for calculation" };
  3843.          stopmsg(0,msg);
  3844.          }
  3845.          break;
  3846.       case STRING1:
  3847.          {
  3848.          static char far msg[]={"String can be a maximum of 10 digits" };
  3849.          stopmsg(0,msg);
  3850.          }
  3851.          break;
  3852.       case STRING2:
  3853.          {
  3854.          static char far msg[]={"Make string of 0's through  's" };
  3855.          msg[27] = k_1 + 48; /* turn into a character value */
  3856.          stopmsg(0,msg);
  3857.          }
  3858.          break;
  3859.       case TABLEK:
  3860.          {
  3861.          static char far msg[]={"Make table of 0's through  's" };
  3862.          msg[26] = k_1 + 48; /* turn into a character value */
  3863.          stopmsg(0,msg);
  3864.          }
  3865.          break;
  3866.       case TYPEKR:
  3867.          {
  3868.          static char far msg[]={"Type must be 21, 31, 41, 22, 23, 24" };
  3869.          stopmsg(0,msg);
  3870.          }
  3871.          break;
  3872.       case RULELENGTH:
  3873.          {
  3874.          static char far msg[]={"Rule must be    digits long" };
  3875.          i = rule_digits / 10;
  3876.          if(i==0)
  3877.             msg[14] = rule_digits + 48;
  3878.          else {
  3879.             msg[13] = i;
  3880.             msg[14] = rule_digits % 10 + 48;
  3881.          }
  3882.          stopmsg(0,msg);
  3883.          }
  3884.          break;
  3885.       case CELLULAR_DONE:
  3886.          break;
  3887.    }   
  3888.    if(cell_array[0] != NULL)
  3889. #ifndef XFRACT
  3890.       cell_array[0] = NULL;
  3891. #else
  3892.       farmemfree((char far *)cell_array[0]);
  3893. #endif
  3894.    if(cell_array[1] != NULL)
  3895. #ifndef XFRACT
  3896.       cell_array[1] = NULL;
  3897. #else
  3898.       farmemfree((char far *)cell_array[1]);
  3899. #endif
  3900. }   
  3901.  
  3902. cellular () {
  3903.    S16 start_row;
  3904.    S16 filled, notfilled;
  3905.    U16 cell_table[32];
  3906.    U16 init_string[16];
  3907.    U16 kr, k;
  3908.    U32 n, lnnmbr;
  3909.    U16 i,j,l;
  3910.    S16 t;
  3911.    S32 randparam;
  3912.  
  3913.    set_Cellular_palette();
  3914.  
  3915.    randparam = param[0];
  3916.    lnnmbr = param[3];
  3917.    kr = param[2];
  3918.    switch (kr) {
  3919.      case 21:
  3920.      case 31:
  3921.      case 41:
  3922.      case 22:
  3923.      case 23:
  3924.      case 24:
  3925.         break;
  3926.      default:
  3927.         abort_cellular(TYPEKR, 0);
  3928.         return -1;
  3929.         break;
  3930.    }
  3931.  
  3932.    r = kr % 10; /* Number of nearest neighbors to sum */
  3933.    k = kr / 10; /* Number of different states, k=3 has states 0,1,2 */
  3934.    k_1 = k - 1; /* Highest state value, k=3 has highest state value of 2 */
  3935.    rule_digits = (r * 2 + 1) * k_1 + 1; /* Number of digits in the rule */
  3936.  
  3937.    if ((!rflag) && randparam == -1)
  3938.        --rseed;
  3939.    if (randparam != 0 && randparam != -1) {
  3940.       n = param[0];
  3941.       t = (S16) log10(n) + 1; /* # of digits in initial string */
  3942.       if (t>10 || t <= 0) {
  3943.          abort_cellular(STRING1, 0);
  3944.          return -1;
  3945.       }
  3946.       for (i=0;i<16;i++)
  3947.          init_string[i] = 0; /* zero the array */
  3948.       t = (S16) ((16 + t)/2);
  3949.       for (i=t;i>0;i--) { /* center initial string in array */
  3950.          init_string[i] = n%10;
  3951.          if (init_string[i]>k_1) {
  3952.             abort_cellular(STRING2, 0);
  3953.             return -1;
  3954.          }
  3955.          n /= 10;
  3956.       }
  3957.    }
  3958.  
  3959.    srand(rseed);
  3960.    if (!rflag) ++rseed;
  3961.  
  3962. /* generate rule table from parameter 1 */
  3963. #ifndef XFRACT
  3964.    n = param[1];
  3965. #else
  3966.    /* gcc can't manage to convert a big double to an unsigned long properly. */
  3967.    if (param[1]>0x7fffffff) {
  3968.        n = (param[1]-0x7fffffff);
  3969.        n += 0x7fffffff;
  3970.    } else {
  3971.        n = param[1];
  3972.    }
  3973. #endif
  3974.    if (n == 0) { /* calculate a random rule */
  3975.       n = rand15()%k;
  3976.       for (i=1;i<rule_digits;i++) {
  3977.          n *= 10;
  3978.          n += rand15()%k;
  3979.       }
  3980.       param[1] = n;
  3981.    }
  3982.    if (rule_digits < (int)log10(n)) { /* leading 0s could make log(n) smaller */
  3983.       abort_cellular(RULELENGTH, 0);
  3984.       return -1;
  3985.    }
  3986.    for (i=0;i<rule_digits;i++) {
  3987.       cell_table[i] = n%10;
  3988.       if (cell_table[i]>k_1) {
  3989.          abort_cellular(TABLEK, 0);
  3990.          return -1;
  3991.       }
  3992.       n /= 10;
  3993.    }
  3994.  
  3995.  
  3996.    start_row = 0;
  3997. #ifndef XFRACT
  3998.   /* two 4096 byte arrays, at present at most 2024 + 1 bytes should be */
  3999.   /* needed in each array (max screen width + 1) */
  4000.    cell_array[0] = (BYTE *)&dstack[0]; /* dstack is in general.asm */
  4001.    cell_array[1] = (BYTE *)&suffix[0]; /* suffix is in general.asm */
  4002. #else
  4003.    cell_array[0] = (BYTE far *)farmemalloc(ixstop+1);
  4004.    cell_array[1] = (BYTE far *)farmemalloc(ixstop+1);
  4005. #endif
  4006.    if (cell_array[0]==NULL || cell_array[1]==NULL) {
  4007.       abort_cellular(BAD_MEM, 0);
  4008.       return(-1);
  4009.    }
  4010.  
  4011. /* nxtscreenflag set to 1 by space bar in fractint.c */
  4012. /* set to 2 by non screen writing loop to restart loop on interrupt */
  4013.  
  4014.    filled = 0;
  4015.    notfilled = 1-filled;
  4016.    if (resuming && !nxtscreenflag) {
  4017.       start_resume();
  4018.       get_resume(sizeof(int),&start_row,0);
  4019.       end_resume();
  4020.       get_line(start_row,0,ixstop,cell_array[filled]);
  4021.    }
  4022.    else if (nxtscreenflag == 1) {
  4023.       start_resume();
  4024.       end_resume();
  4025.       get_line(iystop,0,ixstop,cell_array[filled]);
  4026.       param[3] += iystop + 1;
  4027.       start_row = -1; /* after 1st iteration its = 0 */
  4028.       nxtscreenflag = 0; /* reset flag */
  4029.    }
  4030.    else {
  4031.     if(rflag || randparam==0 || randparam==-1){
  4032.       for (col=0;col<=ixstop;col++) {
  4033.          cell_array[filled][col] = rand15()%k;
  4034.       }
  4035.     } /* end of if random */
  4036.  
  4037.     else {
  4038.       for (col=0;col<=ixstop;col++) { /* Clear from end to end */
  4039.          cell_array[filled][col] = 0;
  4040.       }
  4041.       i = 0;
  4042.       for (col=(ixstop-16)/2;col<(ixstop+16)/2;col++) { /* insert initial */
  4043.          cell_array[filled][col] = init_string[i++];    /* string */
  4044.       }
  4045.     } /* end of if not random */
  4046.     if (lnnmbr != 0)
  4047.       nxtscreenflag = 2;
  4048.      else
  4049.       put_line(start_row,0,ixstop,cell_array[filled]);
  4050.    }
  4051.    start_row++;
  4052.  
  4053. /* This section calculates the starting line when it is not zero */
  4054. /* This section can't be resumed since no screen output is generated */
  4055. /* calculates the (lnnmbr - 1) generation */
  4056.    if (nxtscreenflag == 2) { /* line number != 0 */
  4057.      thinking(1,"Cellular thinking (higher start row takes longer)");
  4058.      for (row = start_row; row < lnnmbr; row++) {
  4059.       if(rflag || randparam==0 || randparam==-1){
  4060.        /* Use a random border */
  4061.        for (i=0;i<=r;i++) {
  4062.          cell_array[notfilled][i]=rand15()%k;
  4063.          cell_array[notfilled][ixstop-i]=rand15()%k;
  4064.        }
  4065.       }
  4066.       else {
  4067.        /* Use a zero border */
  4068.        for (i=0;i<=r;i++) {
  4069.          cell_array[notfilled][i]=0;
  4070.          cell_array[notfilled][ixstop-i]=0;
  4071.        }
  4072.       }
  4073.  
  4074.        for (col=r;col<ixstop-r;col++) {
  4075.          t = 0;
  4076.          for (i=col-r;i<=col+r;i++)
  4077.            t += cell_array[filled][i];
  4078.          if (t>rule_digits || t<0) {
  4079.            thinking(0, NULL);
  4080.            abort_cellular(BAD_T, t);
  4081.            return(-1);
  4082.          }
  4083.          cell_array[notfilled][col] = cell_table[t];
  4084.        }
  4085.        filled = notfilled;
  4086.        notfilled = 1-filled;
  4087.        if (check_key()) {
  4088.           thinking(0, NULL);
  4089.           nxtscreenflag = 2; /* can't resume this, but must start over */
  4090.           abort_cellular(CELLULAR_DONE, 0);
  4091.           alloc_resume(10,1);
  4092.           put_resume(sizeof(int),&row,0);
  4093.           return -1;
  4094.        }
  4095.      }
  4096.    start_row = 0;
  4097.    thinking(0, NULL);
  4098.    nxtscreenflag = 0;
  4099.    }
  4100.  
  4101. /* This section does all the work */
  4102.    for (row = start_row; row <= iystop; row++) {
  4103.  
  4104.       if(rflag || randparam==0 || randparam==-1){
  4105.        /* Use a random border */
  4106.        for (i=0;i<=r;i++) {
  4107.          cell_array[notfilled][i]=rand15()%k;
  4108.          cell_array[notfilled][ixstop-i]=rand15()%k;
  4109.        }
  4110.       }
  4111.       else {
  4112.        /* Use a zero border */
  4113.        for (i=0;i<=r;i++) {
  4114.          cell_array[notfilled][i]=0;
  4115.          cell_array[notfilled][ixstop-i]=0;
  4116.        }
  4117.       }
  4118.  
  4119.        for (col=r;col<ixstop-r;col++) {
  4120.          t = 0;
  4121.          for (i=col-r;i<=col+r;i++)
  4122.            t += cell_array[filled][i];
  4123.          if (t>rule_digits || t<0) {
  4124.            abort_cellular(BAD_T, t);
  4125.            return(-1);
  4126.          }
  4127.          cell_array[notfilled][col] = cell_table[t];
  4128.        }
  4129.        filled = notfilled;
  4130.        notfilled = 1-filled;
  4131.        put_line(row,0,ixstop,cell_array[filled]);
  4132.        if (check_key()) {
  4133.           abort_cellular(CELLULAR_DONE, 0);
  4134.           alloc_resume(10,1);
  4135.           put_resume(sizeof(int),&row,0);
  4136.           return -1;
  4137.        }
  4138.    }
  4139.  
  4140.   abort_cellular(CELLULAR_DONE, 0);
  4141.   return 1;
  4142. }
  4143.  
  4144. int CellularSetup(void)
  4145. {
  4146.    if (!resuming) {
  4147.       nxtscreenflag = 0; /* initialize flag */
  4148.    }
  4149.    timer(0,curfractalspecific->calctype);
  4150.    return(0);
  4151. }
  4152.  
  4153. static void set_Cellular_palette()
  4154. {
  4155.    extern char far *mapdacbox;
  4156.    static Palettetype Red    = {
  4157.       42, 0, 0     };
  4158.    static Palettetype Green  = {
  4159.       10,35,10    };
  4160.    static Palettetype Blue   = {
  4161.       13,12,29    };
  4162.    static Palettetype Yellow = {
  4163.       60,58,18    };
  4164.    int i;
  4165.  
  4166.    if (mapdacbox) return;        /* map= specified */
  4167.  
  4168.    dacbox[0].red  = 0 ;
  4169.    dacbox[0].green= 0 ;
  4170.    dacbox[0].blue = 0 ;
  4171.  
  4172.    dacbox[1].red     = Red.red;
  4173.    dacbox[1].green = Red.green;
  4174.    dacbox[1].blue  = Red.blue;
  4175.  
  4176.    dacbox[2].red   = Green.red;
  4177.    dacbox[2].green = Green.green;
  4178.    dacbox[2].blue  = Green.blue;
  4179.  
  4180.    dacbox[3].red   = Blue.red;
  4181.    dacbox[3].green = Blue.green;
  4182.    dacbox[3].blue  = Blue.blue;
  4183.  
  4184.    dacbox[4].red   = Yellow.red;
  4185.    dacbox[4].green = Yellow.green;
  4186.    dacbox[4].blue  = Yellow.blue;
  4187.  
  4188.    SetTgaColors();
  4189.    spindac(0,1);
  4190. }
  4191.